The Coulomb-Sturmian basis for the nuclear many-body problem 
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Calculations in ab initio no-core configuration interaction (NCCI) approaches, such as the no- 
core shell model (NCSM) or no-core full configuration (NCFC) methods, have conventionally been 
carried out using the harmonic-oscillator many-body basis. However, the rapid falloff (Gaussian 
asymptotics) of the oscillator functions at large radius makes them poorly suited for the description 
of the asymptotic properties of the nuclear wavefunction. We establish the foundations for carrying 
out NCCI calculations with an alternative many-body basis built from Coulomb-Sturmian functions. 
These provide a complete, discrete set of functions with a realistic exponential fallofT. We present 
illustrative NCCI calculations for 6 Li with a Coulomb-Sturmian basis and investigate the center-of- 
mass separation and spurious excitations. 

PACS numbers: 21.60.Cs, 21.10.-k, 27.20.+n, 02.30.Gp 



I. INTRODUCTION 



The combination of powerful theoretical frameworks 
with modern computing capabilities are making possi- 
ble significant advances towards one of the basic goals 
of nuclear theory, namely, an ab inito understanding of 
the nucleus directly as a system of interacting protons 
and neutrons with realistic interactions. Nuclear interac- 
tions motivated by quantum chromodynamics are being 
developed, via effective field theory methods [1, 2], to 
provide an underlying Hamiltonian for the problem. It is 
then necessary to solve the nuclear many-body problem 
for this Hamiltonian, obtaining nuclear eigenstates and 
predictions for observables. In a no-core configuration 
interaction (NCCI) approach, such as the no-core shell 
model (NCSM) [3], the eigenproblem is formulated as 
a matrix diagonalization problem, in which the Hamil- 
tonian matrix is represented with respect to a basis of 
antisymmetrized products of single-particle states. The 
nuclear eigenproblem is then solved for the full A-body 
system of nucleons, i.e., there is no assumption of an 
inert core. 

In practice, NCCI calculations have been based almost 
exclusively on a harmonic oscillator basis. In this arti- 
cle, we consider instead an alternative basis for the NCCI 
approach, built from Coulomb-Sturmian functions [4, 5]. 
These functions have previously been applied to few- 
body problems in atomic [4, 6-8] and hadronic [9-12] 
physics. The Coulomb-Sturmian functions have the dis- 
tinctive property of constituting a complete, discrete set 
of square-integrable functions, while also possessing re- 
alistic exponential asymptotics appropriate to the nu- 
clear problem. In the present work, the foundations for 
carrying out nuclear many-body calculations with the 
Coulomb-Sturmian basis are established. Then, illus- 
trative calculations for the nucleus 6 Li are carried out 



with the NCCI approach in a Coulomb-Sturmian basis. 1 
Many of the considerations addressed here specifically 
in the context of the Coulomb-Sturmian basis are more 
broadly applicable to alternative single-particle bases for 
the nuclear problem. 

Actual NCCI calculations must be carried out in a fi- 
nite, truncated space. Progress in expanding the domain 
of applicability of the method is hampered by a com- 
binatorial scale explosion in the dimension of the prob- 
lem, with increasing size of the included space of single- 
particle states and with the number of nucleons in the 
system. The challenge is to reach a reasonable approxi- 
mation of the converged results which would be achieved 
in the full, untruncated space for the many-body sys- 
tem. The success of the calculation is determined by the 
rate of convergence of calculated observables (energies, 
charge or mass radii, electromagnetic moments and tran- 
sition rates, etc.) with increasing basis size and the abil- 
ity to reliably extrapolate these results for finite spaces 
to the full many-body space [13-15]. Convergence rates 
may be expected to be sensitive to the choice of single- 
particle states from which the NCCI many-body basis is 
constructed, as well as the truncation scheme used for 
the many-body basis. 

Before considering alternative bases, it is worth noting 
that the oscillator functions present significant advan- 
tages as a basis for the nuclear problem, which require 



The Coulomb-Sturmian single-particle states used in the present 
calculations arise as solutions to a general Sturm-Liouville equa- 
tion (Sec. Ill A), rather than a Schrodinger equation or Hartree- 
Fock problem. They consequently do not physically corre- 
spond to "shells" in the conventional sense, i.e., orbitals for 
independent-particle motion in some mean-field potential de- 
scribing the zeroth-order dynamics of the system. Therefore, 
we use the more inclusive term configuration interaction, rather 
than specifically shell model, throughout the present work. 
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further assessment in moving to another basis: 

(1) An exact factorization of center-of-mass and intrin- 
sic wavefunctions is obtained in many-body calculations 
when the oscillator basis is used in conjunction with the 
N max truncation scheme (see Sec. II C), which is based 
on the total number of oscillator quanta. Thus, the os- 
cillator basis with this truncation allows precise removal 
of or correction for spurious center-of-mass contributions 
to the dynamics. 

(2) Matrix elements of the nucleon-nucleon two-body 
interaction are naturally formulated in the relative oscil- 
lator basis, of functions ^ n i(ri— r2) (see Sec. II A). These 
matrix elements can easily be transformed to the two- 
body oscillator basis, of functions ( r i)^n 2 h ( r 2)> by 
the Moshinsky transformation [16]. The simplicity of this 
transformation is lost with any other single-particle ba- 
sis. This is a fundamental concern, since the starting 
point of the many-body calculation is evaluation of the 
two-body matrix elements. (Similar comments apply for 
three-body or higher-body interactions.) 

(3) The oscillator functions constitute a complete dis- 
crete basis for square-integrable functions. Many alter- 
native bases do not provide this convenience. For in- 
stance, the bound state eigenfunctions of the Schrddinger 
equation for finite-depth potentials, such as the Woods- 
Saxon potential, are typically finite in number and in gen- 
eral do not constitute a complete set of square-integrable 
functions, without inclusion of the unbound continuum 
Schrddinger equation solutions as well. 

Nonetheless, there are also strong reasons to consider 
moving beyond the oscillator basis. The classic and long- 
recognized (e.g., Ref. [17]) physical limitation of the oscil- 
lator basis, for application to the nuclear problem, lies in 
the Gaussian falloff (oc e~ ar ) at large distance r, which 
is a consequence of the quadratic confining harmonic os- 
cillator potential. In contrast, for particles bound by a 
finite-range force, the actual asymptotics are exponen- 
tial (cx e~ /3r ). This mismatch in asymptotics, i.e., the 
wavefunction tails, between the expansion basis and the 
physical system imposes a serious handicap on the con- 
vergence of calculations with increasing basis size. The 
problem is especially significant for observables, such as 
the root-mean-square radius or E2 strengths, which are 
sensitive to the large-r properties of the nuclear wave- 
functions. 

To adapt the Coulomb-Sturmian basis to the nuclear 
many-body problem, we must overcome the aforemen- 
tioned technical challenges of moving away from the oscil- 
lator basis. The Coulomb-Sturmian functions, as already 
noted, are complete and offer the convenience of being a 
discrete set. The remaining challenges — transforma- 
tion of matrix elements and center-of-mass factorization 
or spuriosity — are found to be tractable. First, we re- 
view the relevant aspects of the NCCI approach as con- 
ventionally implemented, including the oscillator single- 
particle basis (Sec. II A), the Hamiltonian (Sec. II B), 
and the N max many-body truncation scheme (Sec. II C). 
Then, procedures and results are established for using 



the Coulomb-Sturmian basis for nuclear many-body cal- 
culations. The Coulomb-Sturmian functions are defined 
(Sec. Ill A), practicalities related to the radial length pa- 
rameter are considered (Sec. IIIB), the transformation 
of interaction two-body matrix elements from the oscil- 
lator basis to the Coulomb-Sturmian basis is addressed 
(Sec. IIIC), and it is shown how the two-body matrix 
elements of the relative kinetic energy (and certain other 
operators) can be evaluated via separability (Sec. HID). 
Finally, illustrative NCCI calculations for 6 Li with the 
Coulomb-Sturmian basis (Sec. IV A) are compared with 
oscillator-basis calculations of the same dimensionality. 
The convergence of energies (Sec. IV B) and the root- 
mean-square radius (Sec. IV C) is examined, and issues 
of center-of-mass factorization and spurious states are 
explored in detail (Sec. IV D). Preliminary results were 
reported in Ref. [18]. 



II. BACKGROUND: NO-CORE SHELL MODEL 

A. Harmonic-oscillator basis 

The basis states conventionally used in the NCCI ap- 
proach are antisymmetrized products of single-particle 
harmonic oscillator states. These single-particle states 
are eigenstates of the Hamiltonian 
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where il denotes the oscillator frequency and m N the 
nucleon mass, and r and p are the single-particle coordi- 
nates and momenta. For the spatial part of the solution, 
we have the usual three-dimensional oscillator wavefunc- 
tions 



*»Jm(r) = N nl (r/b) l L l n +1 / 2 [(r/b) 2 ]e-^^^ 2 Y lm (r), 



with normalization factor 
N nl 



1 



2nl 



(Z + n + 1/2)! 
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(2) 



(3) 



where the L" are generalized Laguerre polynomials, the 
Yi m are spherical harmonics, n is the radial quantum 
number, I and m are the orbital angular momentum and 
z-projection, and b is the oscillator length, given by b — 
[^/(raArfi)] 1 / 2 . We use factorial notation [x\ = T(x + 1)] 
uniformly, for both integer and half-integer arguments. 
Letting *„ im (r) = r" 1 R nl (r)Y lm (r), the radial wave- 
function is thus 

Rnl(r) = bN nl {rlb) l + l L l + ll2 [(r/b) 2 ]e-^l^l 2 . (4) 

The functions R n i form an orthonormal set, with 
J °° dr R n 'i(r)R n i(r) — 5 n > n . The full single-particle 
states \nljm), including spatial and spin degrees of free- 
dom, are defined as usual for the nuclear shell model, 
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by coupling the orbital and spin-i angular momenta to 
good total angular momentum j, with its z-projection 
again denoted by m. 

The many-body basis states, for calculations in a space 
of fixed total many-body angular momentum projection 
M (M -scheme basis), are then 

V> = A\n 1 l 1 j 1 m 1 )\n 2 l2j2m 2 ) ■ ■ ■ \n A l A j A m A ) , (5) 

where the operator A represents antisymmetrization, 
over protons and neutrons separately. The many-body 
basis states are thus eigenstates of the Hamiltonian for 
noninteracting particles in a harmonic oscillator poten- 
tial, H n — J^i h?- These states may be classified ac- 
cording to the total number of oscillator quanta iV to t = 
~^2iNi — + h) an d have energy eigenvalue E — 

(N t ot + j)AMl with respect to H n . Thus, truncations 
by iV tot , as considered in the following section, are energy 
truncations under this noninteracting Hamiltonian. 

Since we will later need to consider momentum-space 
wavefunctions, note that these are obtained as the 
Fourier transform 



*„ ;m (k) = (2tt)- 3 / 2 / d 3 re^^ nlm (r). 



(6) 



The radial wavefunction R n i in momentum space is de- 
fined by 



i lm (k) = (-i)'^^(k) 
and is obtained as the Fourier-Bessel transform [19] 



(7) 



Rm(k) = (2/ir) 1 / 2 / drkrji(kr)Rra(r). (8) 
Jo 

For the oscillator, R n i has the same functional form as 
the coordinate-space oscillator wavefunction R n i , with [5] 

Rnl(k) = (-)"^„ ; (6fc)' +1 ^ +1/2 [(&fc) 2 ]e- ( " fc)2/2 , (9) 



where 



Nm = fo 3/2 



2n! 



(Z + n+1/2)! 



1/2 



(10) 



where T is the one-body kinetic energy operator and V 
represents the interaction of the nucleons. Commonly, 
the isoscalar kinetic energy 



n. at ^ 



2m n 



Pi 



(12) 



is used, that is, protons and neutrons are treated cquiv- 
alently as having the average nucleon mass toat, and the 
summation index i runs over all A nucleons. The poten- 
tial V is a Galilean-invariant operator involving two-body 
and possibly higher many-body terms. 

The Hamiltonian (11) has the essential property that 
it may be separated into center-of-mass and intrinsic 
(Galilean-invariant) contributions. The kinetic energy 
operator separates into a term 



2Am 



N 



2Am 



N 



(13) 



representing the center-of-mass kinetic energy and a term 
1 

AAm N 



representing the kinetic energy of relative motion of the 
nucleons, where the prime on the summation Y^!ij indi- 
cates i ^ j. The decomposition of both the r 2 and p 2 
operators into center-of-mass and relative contributions 
is summarized in Appendix A, which also serves to define 
a uniform notation for the present work. The operator 
T rc i depends only upon relative momenta p i — and 
is therefore Galilean invariant. Thus, the full nuclear 
Hamiltonian (11) may be separated as H = T c . m . + H lni 
where 



Hin — T'rel + V 



(15) 



is the Galilean-invariant intrinsic Hamiltonian. As a con- 
sequence of the separability of H , a complete set of eigen- 
states may be found with coordinate-space wavefunctions 
which have the factorized form 



^(r^erj) = Vc.m.(R)V'in,fc(rij;o- i ). 



(16) 



The factor ^ c . m .(R) depends only on the center-of-mass 
coordinate, and the factor i/'in (fij ; o"« ) depends only on 



The Rni form an orthonormal set, with relative coordinates r„ = r; - v, and intrinsic spin de- 



r o °° dkR n >i(k)R n i(k) = 5 n < n . 



B. Hamiltonian 

We now review the properties of the nuclear Hamil- 
tonian which are most relevant to understanding the 
solution method based on the Coulomb-Sturmian ba- 
sis (Sec. Ill) and the results from applying this method 
(Sec. IV). The NCCI approach is based upon a nonrela- 
tivistic nuclear many-body Hamiltonian of the form 



H = T + V, 



(11) 



grees of freedom, indicated schematically here by the ar- 
guments (Ti. For each intrinsic excitation, with wave- 
function V'in, an infinite set of eigenstates sharing this 
same intrinsic structure but different center-of-mass ex- 
citations Vc.m. is obtained. The corresponding energy 
eigenvalue separates into eigenvalues of T c . m . and H m , as 
E = E c , m , + Ei n . 

The "interesting" many-body spectroscopy of the nu- 
cleus resides in the intrinsic wavefunctions ipi n and eigen- 
values Si„, but the "uninteresting" center-of-mass mo- 
tion remains as an unavoidable and potentially obfus- 
cating element of the solution. In principle, the center- 
of-mass motion may be completely eliminated from the 
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problem, by explicitly changing variables to relative co- 
ordinates. However, the nuclear many-body state must 
be antisymmetrized, and this process rapidly becomes in- 
tractable with increasing nucleon number. On the other 
hand, if we instead solve the nuclear eigenproblem in 
a many-body basis constructed from antisymmetrized 
products of single-particle states, antisymmetrization is 
straightforward, but we are consigned to simultaneously 
solving for center-of-mass and intrinsic excitations. 

Before we consider the specifics of formulating the 
eigenproblem with respect to a basis, it is worth consid- 
ering the solutions in the full coordinate space further. 
First, it is convenient to remove the complication of the 
center-of-mass kinetic energy operator, by considering 
the eigenproblem not for the full Hamiltonian H of (11) 
but rather for the intrinsic Hamiltonian H ln of (15). The 
full spectroscopic information of the original problem is 
maintained, since the eigenstates still have wavefunctions 
of the form ^(r^fi) = ?/' c . m .(R)V'in(rij; <Tj), but these 
are now simply associated with eigenvalues E = E ln . 
Thus, for each intrinsic wavefunction V'in, an infinite set 
of eigenstates sharing the same intrinsic structure but 
different center-of-mass excitations ip c . m . is still obtained, 
and these are now strictly degenerate with each other. 

Since T c m . has been eliminated from the Hamilto- 
nian, we are free to consider any complete set of wave- 
functions to span the degenerate space of center-of- 
mass wavefunctions. For instance, suppose plane wave 
solutions ^ c . m .(R) = e * K R arc taken for the cen- 
ter of mass. Then, for each intrinsic excitation ipi n ,k, 
with intrinsic eigenvalue Ek, a continuum of eigen- 
states will be obtained, having wavefunctions erj) = 
e -sK-R^,. n fc ^ r Under the full Hamiltonian H, 
these states form a continuum with E = h 2 K 2 / (2Am^) + 
Ek, but, under H m , these states are infinitely degenerate, 
all with E ~ E k - 

Although they provide the simplest illustration, plane 
wave center-of-mass wavefunctions do not naturally oc- 
cur in our actual solutions to the eigenproblem, which 
are obtained in terms of spatially localized single- 
particle basis wavefunctions. The choice of basis for 
center-of-mass wavefunctions with direct practical sig- 
nificance in oscillator-basis calculations consists instead 
of three-dimensional harmonic oscillator wavefunctions, 
V'c.m.(R) = *„; m (R). The *„ im (R) are eigenfunc- 
tions of the center-of-mass harmonic oscillator Hamilto- 
nian -ff^m. j defined with oscillator frequency f2 and mass 
AniN, i-e., 



point energy of center-of-mass motion with respect to this 
term, so the Lawson term has the form \(H^ m — §M7), 
with A positive, or, more transparently, aN^ m , where 
N^ m — (H^m. — is the number operator as- 

sociated with (17). The Hamiltonian thus becomes 



H = T Ici + V + oJVJ 



Q 

cm. 



(18) 



The factorized eigenstates have coordinate space wave- 
functions ip(Ti; (Ti) = *nZm(R)V'in(rij :; ct;), and the eigen- 
values are now E = Ek + aN c m . Thus, the eigenvalues 
for states with N c m = arc unchanged by the Law- 
son term, still simply the intrinsic energies Ek, while the 
eigenvalues of spurious states, with N c m > 0, are raised 
out of the low-lying spectrum, to an excitation energy of 
at least a. 



C. Many-body iV max truncation 

The factorization of the wavefunction just described 
is possible in the full space of the many-body system. 
However, in practice, diagonalization of the Hamilto- 
nian must be carried out in a finite-dimensional sub- 
space spanned by some truncated basis. In general, one 
cannot expect to be able to construct center-of-mass 
factorized states in such a subspace. The separation 
^(r^crj) = V'c.m. (R)V'in ', i) will be lost, and it will 
not be possible to divide the set of eigenstates into "non- 
spurious" states, consisting of a simple product of a 0s 
center-of-mass wavefunction with a single intrinsic exci- 
tation, and "spurious" states, involving center-of-mass 
excitations. However, there is an important special case 
in which factorization occurs, namely, for a harmonic- 
oscillator many-body basis in the so-called iV max trun- 
cation scheme, which is based on the total number of 
oscillator quanta for the many-body state. This trun- 
cation is commonly used in NCCI calculations. In this 
section, we briefly examine the structure of the -/V max - 
truncated space, both to understand what changes as we 
go to a general single-particle basis and as a prerequisite 
to understanding the spurious state spectrum observed 
for NCCI calculations with the Coulomb-Sturmian basis 
in Sec. IV D. 

Factorization is to be expected if the truncated space 
H for the calculation has a simple product structure, 
before antisymmetrization, 2 



n c.m. ~ ± c. 



+ 



Am N n 2 R 2 



(17) 



n 



U C .m. ® Hi, 



(19) 



The center-of-mass excitation is thus characterized by 
the number iV c m = 2n + I of oscillator quanta. This 
particular choice of center-of-mass wavefunctions is en- 
forced as the eigenstates of the Hamiltonian, if the de- 
generacy of center-of-mass states is broken by introduc- 
ing a Lawson term [20] proportional to H^ m . It is 
both conventional and convenient to subtract the zero- 



2 If the truncated space has the form W c . m . (g) W in , it is in prin- 
ciple possible to choose a basis consisting of product functions 
<Pc.m.,i<Pin,j ■ Since Hj n acts only on intrinsic degrees of free- 
dom, it does not connect basis states involving different (p c .m.,i- 
Therefore, the Hamiltonian matrix with respect to this basis is 
block diagonal, with each block simply consisting of the matrix 
representation of H ln on the basis of intrinsic states (f>in,j- 
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Most simply, if all nucleons are restricted to occupy a 
filled core plus valence orbitals taken from a single ma- 
jor oscillator shell, the many-body space does factorize 
in the form (19), with pure Os motion for the center of 
mass, as shown by Elliott and Skyrme [21]. The essen- 
tial reason is that the total number Atot of harmonic 
oscillator quanta is identical whether calculated as the 
sum of single particle oscillator quanta, A tot = 
or as the sum of center-of-mass and intrinsic quanta 
N tot = A c m + N in [21]. The equivalence may be seen 
from the decomposition of the one-body number oper- 
ator A = {Ml)- 1 [p 2 /{2m N ) + {m N n 2 /2)r 2 - iMl/2] 
into center-of-mass and intrinsic parts (which follows 
from Appendix A). Thus, the space for this situation 
is W — %c.m. ®"^?n> wnere %c. m . i s the one-dimensional 
space containing the Os oscillator function, and Hf n is the 
space of intrinsic functions with no excitations above the 
valence shell. 

The A max truncation scheme is a generalization, which 
likewise yields factorized eigenstates. Consider a space 
spanned by product states subject to the truncation 

N tot = ^2N i <N + N max , (20) 

i 

where Ao is the minimal number of oscillator quanta for 
the given number of protons and neutrons, if all nucle- 
ons occupy the lowest permitted shells. (The Elliott and 
Skyrme space described above is obtained for A max = 0.) 
The A max -truncated space may be decomposed as a di- 
rect sum of product spaces, before antisymmetrization, 3 

+ «c.m. ® " 2 + ' ' ' + W^ST ® (21) 

where is the space of center-of-mass functions with 

exactly A oscillator quanta, and is the space of 
intrinsic functions with A or fewer intrinsic excitation 
quanta above No. Consequently, factorization is main- 
tained, but, in the solution to the many-body problem in 
an A max -truncated space, several approximate copies of 
the intrinsic spectroscopy are obtained, each in a more 
highly-truncated space. The m <8> "Hj^ max block yields 
the "interesting" solutions, or nonspurious states, con- 
sisting of a Os center-of-mass function multiplied by the 
solutions in the least-truncated intrinsic space "Hj^ max . 



3 Since the iVmax-truncated space has the form (21), it is in prin- 
ciple possible to obtain a basis for M^ 111 " 1 consisting of products 

N N . —N 

of the form <f> c ^ m j^; nj - x c m . (The actual basis used in NCCI 
calculations need not be, and generally is not, of this form.) Since 
H ln does not connect basis states involving different center-of- 
mass wavefunctions, the Hamiltonian matrix with respect to this 
basis is block diagonal, with each block, corresponding to a given 
^fjj,*"* , simply consisting of the matrix representation of H ln on 
the basis of intrinsic states for -}{ N ™ 1 ^ X "cm. 



Then the "H^ m <8> 'H 1 ^ max_1 block yields a Op center-of- 
mass function multiplied by the solutions of the intrinsic 
problem in the - H 1 ^ max ~ 1 space, the 'Hc. m . ( 8 ) 'H n ^ max ^ 2 block 
yields Is and Od center-of-mass functions multiplied by 
the solutions of the intrinsic problem in the %^ max_2 
space, etc. In actual calculations, these "uninteresting" 
solutions, or spurious states, may be identified by evalu- 
ating the expectation value (N^ m ) . 

The presence of such spurious states in the low-lying 
calculated spectrum has considerable practical implica- 
tions. Although these states are clearly identifiable, as 
noted, diagonalization of such large matrices as encoun- 
tered in NCCI calculations typically relies upon meth- 
ods such as the Lanczos algorithm [22], which efficiently 
extract a selected set of energy eigenvalues (and corre- 
sponding eigenvectors), namely, those lowest in the en- 
ergy spectrum. With increasing A max , the low-energy 
spectrum would be increasingly cluttered with spurious 
states (as illustrated more concretely in Sec. IV D), lim- 
iting the ability of the Lanczos diagonalization to access 
the low-lying intrinsic excited states. The spurious states 
are therefore, in practice, typically shifted to higher en- 
ergy by inclusion of a Lawson term (Sec. II B) in the 
Hamiltonian, so that they do not interfere with the low- 
lying spectrum obtained by diagonalization. 

As a final practical matter, it is necessary to note that, 
for calculations with parity-conserving nuclear interac- 
tions, the A max truncation of (20) is further restricted 
either to A tot even or to A tot odd. If, e.g., even A tot are 
taken, so - H JVmax is the even-parity space for the nucleus, 
then the %® m ® "H^ max subspace yields only even-parity 
intrinsic excitations, the "Hc. m .(8''H n ^ max_1 subspace yields 
the odd-parity Op center-of-mass function multiplied by 
odd-parity intrinsic excitations, the % 2 m (g)'H 1 ^ max_2 sub- 
space yields the even-parity intrinsic excitations again 
but evaluated in the smaller A max — 2 intrinsic space, 
etc. 



III. THE COULOMB-STURMIAN BASIS 

A. Coulomb-Sturmian functions 

The harmonic oscillator functions have the desirable 
properties, as basis functions for an eigenfunction expan- 
sion, that these form a complete discrete set (of square- 
integrable functions on M 3 ), i.e., without a continuum. 
However, the oscillator functions are obtained from an in- 
finitely bound potential and decay with Gaussian (e _Qr ) 
asymptotics, i.e., they satisfy an undesirable boundary 
condition for problems involving finite binding. Con- 
versely, the Schrddinger equation for the Coulomb po- 
tential yields a set of eigenfunctions which have expo- 
nentially decaying asymptotics (e~^ r ), as desired, but 
which do not form a complete set (of square-integrable 
functions on E 3 ) unless the positive-energy continuum 
Coulomb wavefunctions are included. However, a closely 
related set of functions, the Coulomb-Sturmian func- 
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tions [4-6, 8, 11], can be obtained as the solutions to 
a Sturm-Liouville problem associated with the Coulomb 
potential. These functions retain the exponential asymp- 
totics of the Coulomb problem while also forming, in the 
final form in which we will write them, a complete and 
discrete set of square-integrable functions on W 3 . The 
Coulomb-Sturmian functions thus combine favorable at- 
tributes of both the oscillator and Coulomb functions, 
as an expansion basis for three-dimensional Schrodinger 
problems. 

To begin with, let us recall the Schrodinger equation 
solutions for the Coulomb potential. The functions 



with 



W nlm (r) =N nl 
x L 2 n l+1 

with 

N nl = 



2r 



n + l + 1 
2r 



1 



3/2 



2(n + l + l){n + 2l + l)\ 



satisfy the Schrodinger equation 



~--2E nl )W(r) 



with energy eigenvalue 

E n l = - 



e- r ^ n+l +^Y lm (v), (22) 

1/2 

(23) 
(24) 

(25) 



2(n + / + l) 2 



This is the Schrodinger equation, written in dimension- 
less form (H 2 /m = 1), for the potential V(r) = 1/r. The 
functions W are orthonormal with respect to the stan- 
dard inner product on K 3 , that is, 



d 3 rW:, l , m ,(r)W nlm (r)=5 



(n'V ' m')(nlni) • 



(26) 



Observe that r always appears in the usual Coulomb 
functions divided by a scale n+l + 1, which depends upon 
the quantum numbers n and I. 4 The Coulomb-Sturmian 
functions are obtained by replacing (n+l+1) — > b in (22), 
that is, by carrying out a radial change of variable on each 
function so as to obtain a constant length scale b, yielding 

<S>nlm(r) = N nl (2r/b) l L 2 r !+ 1 (2r/b)e-^ b Y lm (r), (27) 



N„ 



2\3/2 



1/2 



_2(n + / + l)(n + 2/ + l)!_ 

By making the same change of variable in (24), it is seen 
that the functions $ satisfy 

1 

62 



Ctnl- 



$(r) = 0, 



(29) 



with eigenvalue a n i = (n + I + l)/b. They arc thus 
solutions to a Sturm-Liouville eigenproblem, with the 
Coulomb potential as weighting function. 5 The solutions 
$(r) consequently are orthogonal, with respect to the 
same weighting function. In particular, 



, \ 1 

'(r)-$„i. 
r 



1 c 

(30) 

Since (29) is obtained from the Schrodinger equa- 
tion simply by a change of variable, the solutions <fr n im 
may also be considered [4] as a set of solutions to the 
Schrodinger equation. However, by comparison of (24) 
with (29), it is seen then that the scale, or depth, of the 
potential must be taken to vary with each solution, as 
a n i , so the solutions to the problem share a constant en- 
ergy E = — l/(2b 2 ), equal to the ground state energy 
E o of the associated Schrodinger equation, from (25), 
after the substitution (n + I + 1) — > b. 

For use as an expansion basis in quantum mechani- 
cal problems, it is desirable to obtain a set of functions 
which are orthonormal with respect to the standard in- 
tegration metric. This may be accomplished by absorb- 
ing the integration weight 1/r and norm l/[b(n + 1 + 1)] 
appearing in (30) into the Coulomb-Sturmian function 
itself, i.e., multiplying the function $„; m of (27) by 
[b(n + I + Vj/r] 1 / 2 . However, the radial dependence of 
the resulting functions involves a half-integral power of 
r, $ ~ r' -1 / 2 , for r — > 0. In contrast, the harmonic os- 
cillator functions (2) have dependence "J ~ r l for r — > 0. 
We can recover this relation between the r — ¥ asymp- 
totics and the angular momentum by furthermore shift- 
ing I — > 1 + 1/2 in the radial part of the Coulomb- 
Sturmian functions, yielding new functions [5, 26] 

A nlm (r) = N nl {2r/b) l L 2 } +2 {2r/b)e- r ' b Y lm {v), (31) 
where now 



N, 



nl 



-0 



2\3/2 



nl 



(n + 2l + 2)\ 



1/2 



(32) 



4 The combination n + l + 1 is in fact the principal, or energy, 
quantum number, which enters into the energy eigenvalue E nt 
in (25). In comparing with the literature, it should be borne in 
mind that, traditionally, the principal quantum number for the 
Coulomb problem is denoted by n [23], and this notation prop- 
agates to some discussions of the Coulomb-Sturmian functions 
(e.g., Rcfs. [4, 5]). However, consistency with conventional nota- 
tion for the oscillator problem [16] and nuclear shell model [24] 
is strongly desirable in the present context. Hence, we reserve 
the symbol n for the radial quantum number (n = 0, 1, . . .). 



More precisely, the one-dimensional radial equation associated 
with (29), 

' d? ( 1(1 + 1) 1\ 21 . , 

ar z V r A b z J r_ 

obtained by setting ^(r) = r^ 1 ip(r)Yi m (r ), has the form 
of a Sturm-Liouville equation [(d/ dr) p(r) (d / dr) + q(r) + 
Xw(r)]u(r) = [25], with weight function w(r) oc 1/r. 
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Although both $„; m and A„/ m are defined in terms of 
generalized Laguerre polynomials L", the polynomials 
appearing in the $ n ; m have odd a = 21 + 1, while those 
appearing in the A„; m have even a = 21 + 2. The func- 
tions A(r) are orthogonal with respect to the standard 
inner product, i.e., 



rf3 r A*, ; , m ,(r)A„; TO (r) = 5( n 'l'm')(nlm)- 



(33) 



Moreover, they can be shown to form a complete set on 
the space of square-integrable functions on R 3 [5, 27]. 
Letting A n i m (r) = r _1 S n i(r)Yi m (r), the radial wavefunc- 
tion for our Coulomb-Sturmian expansion basis is thus 

S nl (r) = (2/b)- 1 N nl (2r/b) l+1 L 2 J+ 2 (2r/b)e-^ b . (34) 



orthonormal set, 



with 



The S n i form an 
J °° dr S n n(r)S„i(r) = 5„> n . 

The momentum-space representation of the Coulomb- 
Sturmian functions (as for Coulomb functions in general) 
may be evaluated analytically [5, 11]. This property of 
the basis is particularly useful, in the present application, 
for evaluation of matrix elements of the kinetic energy op- 
erators. The momentum-space wavefunction, defined as 
in (6)-(8), is simply expressed in terms of Jacobi polyno- 
mials. If we let A„ im (k) = fc _1 (-i)'S n ;(fc)Y (TO (k), then 



(bk) 



2 + 1 



l]i+2 " 

with normalization factor 



p(Z+3/2,; + l/2) 



{bk) 2 - 1 



(bk) 2 



X nl=2b s,*Wn + 2l + m 



111/2 



The S n i form an 
/ °° dkS n >i(k)S n i(k) = 5„> 



(n + l + ±)\ 



orthonormal set, 



(35) 

(36) 
with 



B. Length parameter 

For any given value of I, the radial wavefunctions 
S n i(r), with n = 0, 1, . . ., constitute a complete and 
orthogonal set on R + , regardless of the choice of length 
scale parameter b in (34). For the full wavefunctions 
A„i m (r) on M 3 , orthogonality of functions with differ- 
ent I quantum numbers is enforced by the Yi m (r) fac- 
tor, regardless of the radial wavefunction. Therefore, the 
choice of length parameter b may be made independently 
for each /-space, and orthogonality of the basis of single- 
particle states on R 3 will still be preserved. 6 



In fact, when spin is introduced in the single-particle basis, a 
distinct value ftjj may be chosen for the length parameter in- 
dependently for each ij'-space, much as different sets of radial 
wavefunctions are obtained for each Ij value in the shell model 
Woods-Saxon basis [24], Different values may also be chosen for 
the proton and neutron spaces. 



The freedom to define distinct 6; , for different values of 
I, appears to be crucial to the present use of a Coulomb- 
Sturmian basis for the nuclear problem. A many-body 
basis built from oscillator wavefunctions has had consid- 
erable past success in providing a reasonable first approx- 
imation to the central portion of the wavefunctions in the 
nuclear problem and also clearly enjoys the advantage of 
complete separability of center-of-mass motion. As we 
introduce the Coulomb-Sturmian basis, we wish to re- 
tain the successes enjoyed by the oscillator basis, to the 
extent possible, while also now providing for exponential 
asymptotics in the tail, or large r, region. 

If b is simply taken independent of /, Coulomb- 
Sturmian radial functions S n i(r) are obtained as shown 
in Fig. 1 (top). For illustration, we use the dimension- 
less value b = 1 for the length parameter. The first four 
radial functions (0 < n < 3) are shown as probability dis- 
tributions \S n i(r)\ 2 , for I = [Fig. 1(a)], I = 2 [Fig. 1(b)], 
and I = 8 [Fig. 1(c)]. These functions may be compared 
with the corresponding radial functions R n i(r) for the 
harmonic oscillator, as shown in Fig. 1 (bottom), again 
taking the dimensionless value b = 1 for the length pa- 
rameter. For the Snii^r), it may be observed that the 
radial probability distribution migrates rapidly to large 
r as I increases. By I = 8, the n — function [Fig. 1(c)] 
shares virtually no overlap with the several lowest-n os- 
cillator functions [Fig. l(i)]. Physically, it is reasonable 
to expect that the success of the oscillator basis in de- 
scribing the central portion of the nuclear wavefunction 
may be lost in such a basis. Convergence of the descrip- 
tion of center-of-mass motion may also be compromised. 
Computationally, there is a purely pragmatic difficulty 
which effectively precludes calculations with such a basis. 
It will be seen in Sec. Ill C that significant overlaps be- 
tween the low-n members of the Coulomb-Sturmian and 
oscillator bases are required, to carry out a change-of- 
basis transformation on the interaction matrix elements 
with reasonable accuracy. 

We therefore seek an alternative prescription for 
which provides a closer alignment of the low-n Coulomb- 
Sturmian basis functions with the harmonic oscillator ba- 
sis functions. A straightforward, though certainly not 
unique, solution is to choose bi so as to align the node 
of the n = 1 Coulomb-Sturmian function, for the given 
value of I, with the node of the n = 1 oscillator func- 
tion, for this same value of I. It is convenient to work in 
this fashion, with nodes rather than, say, maxima, since 
the nodes are given by the zeros of generalized Laguerre 
polynomials [28]. Let x™ s denote the sth zero of the 
generalized Laguerre polynomial L"(x). The condition 
obtained for bi, relative to the oscillator length 6ho> i s 



X 1A 



which yields the simple analytic result 



h 



21 + 3 



(37) 



(38) 
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FIG. 1: Radial basis functions (shown as squared amplitudes), with < n < 3, for I = (left), I — 2 (center), and / = 8 (right): 
(a-c) Coulomb-Sturmian functions S n i(r), with fixed length parameter 6 = 1, (d-f) rescaled Coulomb-Sturmian functions 
S„i(bi;r), with /-dependent length parameter bi given by the prescription (38), and (g-i) harmonic oscillator functions R n i(r), 
with fixed length parameter b = 1 [= 6ho]- The dots mark the location of the node of the n = 1 function in each panel, and 
the connector lines highlight the shift in this node between the top and middle rows. 



Thus, e.g., b a /b HO « 0.8165, 6i/6 H o « 0.6325, and 
i)2/&HO ~ 0.5345. The nodes under consideration are 
marked by dots in Fig. 1. Selecting bi/bno according 
to (38) yields radially rescaled Coulomb-Sturmian func- 
tions as in Fig. 1 (middle). These functions are seen to 
provide a much closer match to the oscillator functions 
of Fig. 1 (bottom) in the small-r central region, than do 
the unsealed functions of Fig. 1 (top), while still retain- 
ing greater support than the oscillator functions in the 
large-r tail region. 

The optimal approach to choosing the bi may be ex- 
pected to depend upon the problem at hand — nucleus, 
interaction, states of interest, obscrvablcs of interest, and 
many-body truncation scheme in use — and warrants 
thorough investigation. The prescription (38) would ap- 
pear to be a reasonable starting point and is therefore 
used in the example NCCI calculations of Sec. IV. How- 



ever, it remains to be determined what prescription for 
bi might ultimately yield the most rapid convergence in 
the many-body problem. Under some circumstances, it 
may even be appropriate to choose the bi separately for 
the proton and neutron spaces, for instance, for neutron 
halo nuclei. 



C. Transformation of matrix elements 

For the many-body problem, we now consider a ba- 
sis built up from the Coulomb-Sturmian functions A„; m , 
combined with spin to give nlj states as usual. The an- 
gular and spin dependence is thus the same as for the 
harmonic oscillator single-particle states, but with the 
harmonic oscillator radial wavefunctions R n i replaced by 
the S n i. Many-body basis states may be built as anti- 
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symmetrized products of these single-particle states ex- 
actly as before, i.e., according to (5). For the many- 
body calculation, it is necessary for one to evaluate the 
matrix elements of the Hamiltonian with respect to the 
many-body basis states. However, the specific choice of 
single-particle basis enters into the problem only through 
the two-body matrix elements of this Hamiltonian, if the 
interaction is limited to two-body contributions, or three- 
body matrix elements if a three-body interaction is con- 
sidered, etc. Here we consider specifically two-body inter- 
actions and matrix elements, but the discussion readily 
generalizes to higher-body interactions. 

If the two-body matrix elements of an interaction are 
known with respect to the oscillator basis, matrix ele- 
ments with respect to the Coulomb-Sturmian basis may 
then be obtained by a straightforward sum over two-body 
states. Strong practical considerations suggest first gen- 
erating the nuclear interaction two-body matrix elements 
in the oscillator representation. By Galilean invariance, 
the interaction itself is a function only of the relative 
r 2 — ri degree of freedom and the intrinsic spins. Con- 
ventionally, for NCCI calculations, the two-body inter- 
action is first represented via its matrix elements in a 
basis of harmonic oscillator states in the relative spatial 
degree of freedom, coupled to the spins, i.e., \nl;SJ). 
The transformation from a relative oscillator basis to 
a single-particle oscillator basis, i.e., to product states 
\n a laja,nbhjb', J) for the two-particle system, can then 
be carried out through the well-developed framework of 
the Moshinsky transformation [16]. Such a convenient 
means of transformation is not, in general, available for 
other bases. 7 Therefore, only after this transformation 
to single-particle degrees of freedom do we carry out the 
transformation to the Coulomb-Sturmian basis. 

For purposes of discussing the change of basis, let us 
label single-particle orbitals for the oscillator basis by 
unbarred symbols a = (n a l a j a ), b = {riblbjb), etc., and 
those for the Coulomb-Sturmian basis by barred sym- 
bols a = {n a l a j a ), b — (fibhjb), etc- Then the two-body 
matrix elements in the oscillator basis are of the form 
(cd; J\V\ab; J), and we wish to obtain transformed ma- 
trix elements (cd; J|V|ao; J). The basic ingredient is the 
transformation of single-particle states, 

1 6) = ^2(a\a) \a). (39) 

a 

The angular functions Yj TO and the coupling with spin 
to yield j are identical for both bases, so (a\a) = 
(■^nJ a \^n a l a )^(i a ] a )(i a j a )i an d the sum over orbitals a in 
fact only involves a sum over radial quantum numbers 
n a . In writing out the overlap (R na l a \Sn a i a ), it is worth- 
while to explicitly indicate the different choices of length 



We note, however, that the weakly-convergent two-center expan- 
sion methods of Rcf. [5] might provide a viable approach for 
carrying out such a transformation. 




1 2 3 4 5 6 7 



r 

FIG. 2: Integrand R n l(brio', r)Sni{bi; r) for the overlap inte- 
gral (40), taken for a representative case (I = 0, n = 5, and 
n = 5). For this plot, bi/buo is given by the prescription (38), 
and &ho is taken to be unity. 

parameter appearing in R n i(r) and S n i(r), for which we 
adopt the notations R n i(b;r) and S n i(b;r). Then, the 
overlap is given by the radial integral 8 

poo 

(Rni\Sni)= drR n i{buo\r)Sni{bv,r). (40) 
Jo 

Equivalently, the overlaps may be evaluated in momen- 
tum space, as 

/>oo 

(Rni\Sm)= dkR n i{buo;k)S nl {bi;k). (41) 
Jo 

When larger values for the radial quantum numbers 
n or n arc considered, the integrand appearing in the 
overlap integral (40) or (41) is highly oscillatory, as illus- 
trated in Fig. 2 — as is to be expected for overlap inte- 
grals of functions with large numbers of nodes. There- 
fore, care must be taken in evaluating the overlap integral 
through numerical quadrature. Conventional quadrature 
formulas are found to be slowly-converging and unreli- 
able. However, the zeros of the integrand are easily deter- 
mined, from the zeros of the generalized Laguerre polyno- 
mials or Jacobi polynomials, in terms of which the radial 
functions arc defined, as summarized in Appendix B. In- 
tegration can then be carried out in a numerically robust 



The oscillator wavefunctions as defined in (2) arc positive at 
the origin, i.e., as r — s- 0. The Coulomb-Sturmian functions as 
defined in (31) have this property as well. It should be noted 
that a conventional phase factor (— ) n may be included in the 
definition of ^ n i m , so that the functions arc instead positive 
at infinity, i.e., as r — 5- oo. If so, this sign must be accounted 
for in evaluating the transformation bracket (40) for the change 
of basis. Alternatively, the phase convention for the Coulomb- 
Sturmian basis may be adjusted analogously. 
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fashion if the full integration range [0,oo) is first broken 
into intervals between successive zeros. Within each in- 
terval, the integrand is well-behaved, and conventional 
numerical quadrature can be carried out reliably. The 
results may then be summed to give the full integral. It 
is found that a 32-point Gauss-Legendre quadrature on 
each interval suffices for present purposes, yielding nu- 
merical errors of < 10 -8 (and generally much better) for 
calculations involving radial wavefunctions with n < 20. 
Integration in the tail region, between the last zero of the 
integrand and infinity, requires special treatment, since 
Gauss-Legendre quadrature is only defined on finite inter- 
vals. One can map the tail region onto a finite interval by 
a suitable transformation of integration variable. Alter- 
natively, and most simply, the integration may be trun- 
cated at a sufficiently large cutoff r max , e.g., r max /b w 50 
is found to suffice in the present calculations. 

For proton-neutron matrix elements, the two-body 
states transform as 

\ab; J) pn = ^2(a\a)(b\b) \ab; J) pn , (42) 

ab 

and the matrix elements consequently transform as 

(cd; J\V\ab; J) pn 

= J2(a\a){b\b){c\c){d\d) (cd; J\V\ab; J) pn . (43) 

abed 

As noted above for (39), the sums over orbitals a, b, c, 
and d need only traverse the radial quantum numbers 



n a , ni,, n c , and nj, preserving the same angular quantum 
numbers. 

For proton-proton or neutron-neutron matrix ele- 
ments, normalization considerations related to antisym- 
metrization must be taken into account in carrying out 
the transformation. Since different normalization con- 
ventions arise in the description of two-particle states, 
the present conventions are briefly summarized in Ap- 
pendix C. It is easiest to state the transformation rule 
if the two-body matrix elements are defined in terms of 
the antisymmetrized (AS) two-particle states \ab; JM) as 
of (C3), which are properly normalized except in the case 
in which both particles occupy the same orbital. Then, it 
maybe be seen [e.g., by carrying out a change of basis on 
the creation operators in (C3) [29]] that we simply have 

\db; J) AS = Y,(a\a)(b\b) \ab; J) as- (44) 

ab 

Consequently, for the two-body matrix elements, 

(cd; J\V\ab; J) as 

= J2(a\a)(b\b)(c\c)(d\d) (cd; J\V\ab; J) AS . (45) 

abed 

The corresponding expression for the transformation in 
terms of the strictly normalized antisymmetrized (NAS) 
states \ab; JM)nas of (C5) is less transparent, since the 
case of identical orbitals must be treated specially within 
the sum, giving 



(cd;J\V\ab;J) NAS = (1 + 4b)" 1/2 (l + *o*)" 1/2 E ^ + *a6) 1/2 (l + S cd ) 1 / 2 (a\a)(b\b)(c\c)(d\d) (cd; J\V\ab; J) NAS . 

abed 

(46) 

I 



It is trivial to convert between AS and NAS matrix ele- 
ments, and thus to use either relation (45) or (46), but 
it is important to note the distinction. 

For actual calculation of the transformed matrix ele- 
ments, the infinite sums over orbitals appearing in the 
transformation rule (43) and (45) [or (46)] must be trun- 
cated, limited in practice by the available set of oscillator- 
basis matrix elements. If a shell-based cutoff, i.e., by 
number of oscillator quanta, is applied to the single- 
particle space, then N < N cut for the single-particle 
states, and the sum J2 a bcd appearing in (43) and (45) is 

truncated to X^cd^'^'^"^ 11 * ■ For example, the set of 
oscillator basis two-body matrix elements required for a 
transformation with cutoff N cut = 13 (14 shells) consists 
of 9.2 x 10 7 proton-neutron two-body matrix elements 
and 2.3 x 10 proton-proton or neutron-neutron matrix 



elements. 9 The summations (43) or (45) only involve 
matrix elements sharing the same angular momentum J, 
parity P, and isospin projection T z (pn, pp, or nn), and 
thus in practice the transformation may be carried out 
separately for each sector of matrix elements, character- 
ized by these quantum numbers. After transformation, 
substantially fewer matrix elements are required for an 
A max -truncated many-body calculation in the same num- 
ber of shells, e.g., for p-shell nuclei, an N max = 12 calcu- 



These are the possible nonzero two-body matrix elements (Ap- 
pendix C), with single-particle states taken from 14 shells, for 
an interaction which is parity-conserving but with no further as- 
sumptions about isospin (or charge) symmetry. Actual nuclcon- 
nucleon interactions may in fact contain fewer independent ma- 
trix elements. 
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lation involves 14 shells but only 4.1 x 10 6 proton-neutron 
two-body matrix elements and 1.0 x 10 6 proton-proton 
or neutron- neutron matrix elements, due to the further 
restriction on N to t- 

The accuracy of the resulting two-body matrix ele- 
ments obtained for the Coulomb-Sturmian basis depends 
on the inclusion of an adequate number of oscillator 
shells. The effect of truncation may in general be ex- 
pected to vary depending on the two-body operator un- 
der consideration. In practice, the adequacy of the trans- 
formation may be judged by the sensitivity of the fi- 
nal many-body calculation to N cu t. Calculations with 
A cut = 9 (10 shells), 7V cut = 11 (12 shells), and A cut = 13 
(14 shells) are considered in Sec. IV. 

Since the change of basis (39) represents a transfor- 
mation of radial wavefunctions, the underlying approx- 
imation in applying a cutoff is that we arc effectively 
representing the Coulomb-Sturmian radial functions in 
terms of a truncated set of oscillator radial functions, as 

N<N cM 

Sm(bi;r)= ]T (R n i\S n i) R n i(bno;r), (47) 

n 

with N = 2n + I, so n < (N cut — The decomposi- 

tion of Coulomb-Sturmian functions in terms of oscillator 
functions, shown as squared amplitudes (probabilities), is 
given in Fig. 3. Results are shown for the functions pre- 
viously plotted in Fig.l(d-f), that is, with < n < 3 
and for I = 0, 2, and 8, with the length scales of the 
functions determined according to the prescription (38). 
While the first two or three CS functions for each value of 
I are easily expanded in the oscillator basis, the required 
number of shells is seen to grow rapidly for higher radial 
quantum numbers. The degree to which the Coulomb- 
Sturmian radial function is successfully expanded in a 
truncated set of oscillator radial functions is seen from 
the dashed curves in Fig. 3, which indicate the accumu- 
lated probability P n = I]„,<„(-Rn';|S'nz) 2 - The set of 
oscillator radial functions retained in the most generous 
truncation used in Sec. IV, -/V cut = 13, can be seen from 
the vertical dotted line in each panel of Fig. 3. 



D. Evaluation of two-body matrix element for 
separable radial and kinetic operators 

If the two-body matrix elements of the entire Hamilto- 
nian are first evaluated in the oscillator basis then trans- 
formed to the Coulomb-Sturmian basis, according to the 
procedure of Sec. IIIC, it is found (Sec. IV) that the ki- 
netic energy term requires an unacceptably large number 
of oscillator shells for its expansion. That is, the A cut - 
dependence of the transformed relative kinetic energy, 
rather than of the transformed nucleon-nucleon interac- 
tion, dominates the cutoff dependence of the many-body 
calculations. 

In this section, we therefore instead consider a scheme 
which permits the two-body matrix elements of the 



center-of-mass and relative components of the r 2 and p 2 
operators — R 2 , r 2 cl , P 2 , and p 2 cl — to be evaluated 
directly in the Coulomb-Sturmian basis. The approach 
makes use of separability, together with the explicitly 
known form (35) of the Coulomb-Sturmian radial wave- 
function in momentum space. The operators R 2 , r 2 cl , P 2 , 
and p 2 el all appear in the NCCI problem. Specifically, 
p 2 el appears through the relative kinetic energy operator, 
r rei through the root-mean-square (RMS) radius observ- 
able, and R 2 and P 2 through the center-of-mass oscil- 
lator Hamiltonian appearing in the Lawson term. The 
definitions of and relations among these operators are 
summarized for reference in Appendix A. 

Each of the operators R 2 , rf el , P 2 , and p 2 cl may be 
decomposed into one-body terms and separable two-body 
terms. In the following, we let p = hk and work with 
K 2 and k 2 cl instead of P 2 and p 2 cl . Then we have (see 
Appendix A): 

i ij 

A 2 r 2 cl = (A-l)^-Y!^-^ 

i ij 
i ij 



(48) 



k 2 



rel 



I [ i terms involving ^ i arc manifestly one-body op- 
erators, and those involving J2ij are manifestly two- 
body operators. The important property of these expres- 
sions (48) for the present approach is that the two-body 
term in each case — Yl'ij r i ' r j or Yl'ijPi ' Pj — nas the 
separable form X^Ti-Tj, where is a spherical tensor 
(in the present case, rank-1 or vector) operator acting on 
particle k only. The procedure for calculating two-body 
matrix elements therefore reduces to the evaluation of ra- 
dial integrals (either in coordinate space or momentum 
space, for Yi or k i; respectively), which are then com- 
bined using standard angular momentum coupling and 
recoupling results. 

First, let us consider the matrix elements of the one- 
body terms J2i r i an d J2i ^i appearing in (48). For the 
Coulomb-Sturmian basis, the one-body matrix elements 
of r 2 and k 2 are 



(b\r 2 \a) =6 h iJ jt 



and 



(b\k 2 \a) = Si b i a 5j 



•f 

Jo 

r 

Jo 



drS nb i b (bi b ;r)r S nJa {b la ;r) 

(49) 

dk S nb i b {b h ; k) k 2 S nJa (b u ; k) . 

(50) 

The radial integrals appearing in these expressions may 
be evaluated by numerical quadrature. Since the inte- 
grands are highly oscillatory, the comments and meth- 
ods of Sec. Ill C apply to this integration. The integrals 
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FIG. 3: Probability decomposition of the Coulomb-Sturmian radial functions S n i(bf,r) with respect to the basis of harmonic 
oscillator radial functions R„i(buo\ r). Results are shown for Coulomb-Sturmian functions with < n < 3 (top to bottom) and 
for I — (left), I = 2 (center), and I = 8 (right), with bi/buo given by the node- matching prescription (38). The histogram bars 
indicate squared amplitudes (R„i (&ho; r)\Sm(bi; r)) 2 with respect to individual oscillator basis functions. The dashed curve 
indicates accumulated probability, i.e., for all oscillator basis functions of lesser or equal n. The vertical dotted line indicates 
the truncation of the radial basis in effect if the oscillator functions are limited to 14 major shells (JV cut = 13), as in the 
least-truncated calculations of Sec. IV. 



are again evaluated piecewise between zeros of the inte- 
grands, through Gauss-Legendre quadrature. 

Although ^2 1 rf and kf are one-body operators, 
they are being considered here as contributions to the 
two-body operators R 2 , r 2 cl , P 2 , and p 2 el , through (48), 
for which two-body matrix elements are therefore required 
as input to the many-body calculation. The appropriate 
two-body matrix elements are readily obtained from the 
one-body matrix elements (b\r 2 \a) and (b\k 2 \a) consid- 
ered in (49) and (50). In general, corresponding to any 
one-body operator U = itj, we may define a two-body 
operator Vu via Vu — ^Yl'ij v ij' where Vij = Ui + Uj. By 
comparing the sums appearing in the definitions of U and 
Vu, it may be seen that these operators are identical, ex- 
cept for an A-dependent normalization. Specifically, the 



operators are related by 

U =A^l Vu ' (51) 

when acting on the many-body states of an A-particle 
system. 

We therefore consider two-body matrix elements of Vu- 
For the proton-neutron matrix elements, 

(cd; J\Vu\ab; J) pn = (c\U\a)5 db + (d\U\b)5 ca . (52) 

For the proton-proton or neutron-neutron matrix ele- 
ments, the antisymmetrized matrix element may be eval- 
uated by first reexpressing it in terms of unsymmetrized 
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matrix elements, as 

(cd; J\Vu\ab; J) as = {cd; J\vi 2 \ab; J) 

-(-y-i«-io(cd;J\v 12 \ba;J), (53) 

with denned above. It follows that 

(cd; J\Vu\ab; J) as = (c\U\a)S db + (d\U\b)S ca 

- (-) J - j °- jb (c\U\b)*da - (-) J -^ b (d\U\a)S cb . (54) 

Thus, the matrix elements of interest for the one-body 
terms appearing in (48) are obtained by setting U = r 2 
or k 2 and using one-body matrix elements (49) or (50), 
respectively, in (52) and (54). 

Now let us consider the matrix elements of the two- 
body terms lY^'ij r i' r j an d appearing in (48). 
We include a factor of 1/2 in these expressions to bring 
them into the standard form for two-body operators, 
namely, V — ^2ijVij, with = v^. The operator de- 
fined by the sum, in either case, is of the separable form 
Vti-t 2 = sSijTj -Tj, where Tj is a vector operator act- 
ing on particle i. Since the summand is a spherical tensor 
product of operators acting on two different subsystems 
(namely, particles i and j), it is possible to evaluate the 
matrix elements by Racah's reduction formula [30]. For 
the proton-neutron matrix elements, 



(cd;J\T 1 -T 2 \ab;J) pn 



jc jd J 
jb ja 



{}<c||T||a)<d||T||6). (55) 



For the proton-proton and neutron-neutron matrix ele- 
ments, it is important to note that Racah's reduction 
formula applies to matrix elements between ordinary, 
unsymmetrized product states of distinguishable subsys- 
tems. Thus, the two-body matrix element between an- 
tisymmetrized states of two like nuclcons must first be 
expanded by (C3) in terms of unsymmetrized matrix el- 
ements, as 

(cd;J\V Tl .T 2 \ab; J) As = (cd; J\T 1 ■ T 2 \ab; J) 

- (-) J -^- 3b (cd;J\T 1 -T 2 \ba;J). (56) 

Then, each of the two terms may be evaluated separately 
through Racah's reduction formula, much as in (55), giv- 
ing 

(cd; J|Ti • T 2 |a6; J) 

= (-) jMJ {%f a i}(c\\T\\a)(d\\T\\b) (57) 



for the first term, and similarly with b o a for the second 
term. 

The one-body reduced matrix elements (6||T||a) ap- 
pearing in (55) or (57) are expressed in terms of radial in- 
tegrals, using the general relation x m — y/4Tr/3xY lm (x) 
for the spherical components of a coordinate vector x in 
terms of Y\ [24] , as 



<6||r||a> 



47T 

T 



1/2 



/>oo 

/ dr S rlb i b (bi b ;r)r S nJa (b la ;r) 
Ja 

X (hjbWY^laja) (58) 



and 



(6||k||a) = (-) 



= (\(h-l a -l)/2 



47T 



1/2 



j dkS nh i b {b h ;k)kS na i a (b la ;k) 
Jo 

X (l b jb\\Yl\\laja). (59) 



Numerical evaluation of these radial integrals is again 
subject to the considerations for oscillatory integrands 
discussed in Sec. IIIC. The angular factor appearing in 
(58) and (59) is given by [24] 



(hjb\\Yl\\laja) 



_3_ 

An 



1/2 



(-)" +1 (j4lO|j 6 ±)7r(U/ 6 ), 

(60) 

where7r(ZiZ 2 ---) = |[l+(-)' 1+i2+ '"]- The factor tt(Z 1Z 6 ) 
enforces the parity selection rule for Y\, namely, lb — l a 
odd. Since the angular momentum triangle inequal- 
ity also applies, the radial matrix elements (6||r||a) or 
(6||k||a) need only be evaluated for pairs of orbitals for 
which 4 = l a ±l. The phase factor (-)(h-U-i)/2 j n (59) 
arises from the phase factor (— i) in the definition (7) of 
the momentum-space radial wavefunction, after simplifi- 
cations are caried out making use of the constraints on 
Z-values imposed by the angular factor (60). 

To summarize, the two-body matrix elements of R 2 , 
r 2 el , K 2 , or k 2 cl are evaluated by calculating the one 
body contributions according to (52) or (54) and com- 
bining these with the matrix elements of the two-body 
contribution, calculated according to (55) or (56), via 
the operator relations (48). Collecting the various con- 
tributions and normalization factors, we have 
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(cd; J\A 2 R 2 \ab; J) = -(cd; J\V r 2\ab; J) + 2(cd; J\V ri . T2 \ab; J) 

I\ _L 

(cd; J\A 2 r 2 cl \ab; J) = (cd; J\V r 2\ab; J) - 2(cd; J\V ri . r2 \ab; J) 



(cd;J\K 2 \ab; J) 



1 



A-l 



(cd; J\V k 2 \ab; J) + 2(cd; J\V kl . k2 \ab; J) 



(cd; J\k 2 cl \ab; J) = (cd; J|Vfc2 \ab; J) — 2(cd; J|Vki-k 2 \ a b; J). 

I 



(61) 



Further practical aspects of evaluating these matrix ele- 
ments are considered in Appendix D. 

Although the separable method described here for eval- 



uating two-body matrix elements of R 



,, K 2 , and 



k 2 el has been presented in the context of the Coulomb- 
Sturmian basis, this approach is applicable to a general 
radial basis, so long as both the coordinate-space and 
momentum-space radial wavefunction can be accurately 
evaluated and integrated. The only basis dependence 
lies in evaluating the radial integrals (49), (50), (58), 
and (59). For instance, the separable method can be 
used with the oscillator basis, applied to the radial func- 
tions R n i(r) of (4) and R n i(k) of (9), in lieu of Moshinsky 
transformation. 10 



IV. COULOMB-STURMIAN CALCULATIONS 
FOR 6 Li 

A. Overview 

As a basic illustration of the use of the Coulomb- 
Sturmian basis for NCCI calculations, we consider the 
nucleus 6 Li. The code MFDn [31-33] is used for the 
many-body calculations, taking as its input Hamiltonian 
two-body matrix elements obtained according to the pro- 
cedures developed in Sees. IIIC and HID. Calculations 
are carried out with respect to a proton-neutron M- 
scheme basis. 

The question arises as to how to truncate a many- 
body basis built from Coulomb-Sturmian functions. For 
the present calculation, we formally carry over the A^ max 
truncation scheme to the Coulomb-Sturmian basis. That 
is, for each Coulomb-Sturmian single-particle state, we 
define N = 2n + I. Then, as for the oscillator basis, we 
label the many-body states by the sum iVtot = 
and apply the -/V max truncation as defined in (20). Since 
n is now the radial quantum number for the Coulomb- 
Sturmian functions, the label N no longer has any direct 



10 In fact, the separable method has been used to evaluate the 
matrix elements of the T ro i , N c m , and r^ el operators for the 
oscillator-basis NCCI calculations shown in Sec. IV. Compari- 
son against the results obtained with existing Moshinsky-based 
oscillator-basis calculations provides a vital means of validating 
the present computational framework for general bases. 



significance in terms of oscillator quanta. Furthermore, 
when applied to the Coulomb-Sturmian basis, the A^ max 
truncation does not imply the exact separation properties 
described in Sec. II C, nor can it any longer be interpreted 
as an "energy" truncation, with respect to some nonin- 
teracting Hamiltonian. Nonetheless, as one of many con- 
ceivable truncation schemes, the A^ max scheme provides 
a reasonable starting point for further exploration, and 
it is particularly convenient for use with existing NCCI 
many-body codes. Furthermore, using an N max trunca- 
tion facilitates comparison of convergence rates obtained 
using the oscillator and Coulomb-Sturmian bases, since 
the dimensions of the many-body spaces are then the 
same in both cases. 

The result for any given observable has a twofold de- 
pendence on the basis used: on the truncation and on 
the length parameter. In the existing literature on the 
NCCI approach with the oscillator basis, the oscillator 
length b for the basis is commonly not stated directly, 
but rather the oscillator energy Ml is given, in terms of 
which we recall b = [h/ ' (tun^I)] 1 ! 2 . For consistency, we 
therefore adopt the same convention for the Coulomb- 
Sturmian basis. However, it must be borne in mind that 
the Ml value quoted for the Coulomb-Sturmian basis is 
simply the Ml of the reference oscillator length &ho , from 
which the actual ^-dependent length parameters 6; are 
derived by the node-matching prescription of Sec. IIIB. 
It therefore has no direct significance as an energy scale 
for the problem. When comparing calculations in the 
harmonic oscillator basis and in the Coulomb-Sturmian 
basis, the relationship of Ml values between the two cal- 
culations should therefore also not be viewed as one of 
strict physical equivalence, e.g., it is not necessarily most 
appropriate to compare an Ml = 20 McV oscillator ba- 
sis calculation with an Ml = 20 MeV Coulomb-Sturmian 
basis calculation. Rather, a set of calculations for each 
basis, spanning a range of Ml values, should be consid- 
ered, and best convergence may be obtained for different 
Ml values in each of the two bases. However, for either 
basis, the same proportionality b cx (Ml)^ 1 / 2 holds, e.g., 
a doubling in Ml corresponds to a factor of \[2 contrac- 
tion of the length scale. 

The present 6 Li calculations are carried out for the 
JISP16 interaction [34], which is a two-body interac- 
tion derived from neutron-proton scattering data and 
adjusted via a phase-shift equivalent transformation to 
describe light nuclei without explicit three-body interac- 
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tions. All calculations shown here are for the positive- 
parity space, spanned by states with even values N tot — 
N , N + 2, N + iV max . Although isospin is 

not strictly conserved by the Hamiltonian, due to the 
Coulomb interaction, the isospin T is essentially a good 
quantum number for the states in the present calcula- 
tions. Therefore, for simplicity, we restrict attention to 
the T = spectrum. Calculations are carried out in 
several truncated spaces with N max < 10, to provide an 
initial investigation into convergence. 

The nucleus 6 Li provides a useful case for benchmark- 
ing, since calculations with comparatively large values 
of A^max are feasible with the most-powerful presently- 
available computational resources, and detailed extrapo- 
lation studies have recently been carried using the con- 
ventional oscillator basis in such large spaces, specifically, 
JV max < 16, with the same interaction as used here [35]. 
These results provide estimates for the true values of ob- 
servables, against which the present Coulomb-Sturmian 
calculations in smaller spaces can be compared. 



B. Energies 

We begin by comparing the ground state energy ob- 
tained in NCCI calculations with the conventional oscilla- 
tor basis and with the Coulomb-Sturmian basis. The cal- 
culated energies of the 1 + ground state of 6 Li are shown 
for the oscillator basis in Fig. 4(a) and for the Coulomb- 
Sturmian basis in Fig. 4(b). In each case, the calculations 
span a range of Ml values from lOMeV to 40MeV and 
are carried out for iV max = 4, 6, 8, and 10. These iV max 
values correspond to the highest to lowest curves, respec- 
tively, in the figure. The bare Hamiltonian has been used, 
without renormalization to the finite space, so the vari- 
ational principle is in effect, and energies (for the low- 
est state with each set of conserved quantum numbers) 
approach the full-space value monotonically from above 
with increasing N max . 

The goal is not for any single NCCI calculation to ac- 
tually reach a converged value, but rather to obtain the 
most reliable extrapolation from a series of NCCI calcu- 
lations, to the converged value which would be obtained 
in the full, untruncated space for the many-body prob- 
lem [13-15]. It is thus first necessary to examine the de- 
pendence of the results on the basis parameters Ml and 
iV max as just described. Extrapolation schemes are still 
largely empirical in their justification, and different pre- 
scriptions, varying in their details, might be used. How- 
ever, for energies, at least, the basic procedure explored 
in, e.g., Refs. [13, 14, 36], consists of an exponential ex- 
trapolation. The no-core full configuration (NCFC) ap- 
proach [14], in particular, is based on exponential extrap- 
olations of results of calculations obtained with an un- 
renormalized interaction appropriate to the infinite, un- 
truncated space, so that energies approach the full-space 
values monotonically, as noted above. One first finds the 
variational minimum with respect to Ml, for the highest 



available A^x-truncated space. Then one extrapolates 
with respect to iV max , at this Ml, to the full-space result 
(A^max — > oo) by assuming an exponential approach to 
the asymptotic value E^, 

E(N m ^) = E 00 +ae- cN —, (62) 

where E^, a, and c are taken as parameters. 

As the baseline for comparison, the calculations of the 
ground state energy with the oscillator basis are shown 
in Fig. 4(a). The variational minimum with respect to 
Ml occurs at ~ 20McV, for N max = 10, moving gradu- 
ally lower with increasing N max . For each value of Ml at 
which calculations have been carried out, an exponential 
extrapolation of the N mliX — 4-10 calculations is shown 
(indicated by a cross). The best estimate of the ground 
state energy from Ref. [35], E = -31.49(3) MeV, is in- 
dicated by the dashed horizontal line. The extrapolated 
values pass through this estimate at Ml rts 20 MeV, that 
is, roughly the location of the variational minimum. 

The calculations of the ground state energy with the 
Coulomb-Sturmian basis are shown in Fig. 4(b). The 
variational minimum with respect to Ml occurs at ~ 
30 MeV, for A^x = 10, and moves higher with increas- 
ing A^max- Notice that at each N max the variational min- 
imum energy obtained with the Coulomb-Sturmian basis 
is substantially higher than that obtained with the os- 
cillator basis (by ~ 2 MeV for N max = 10). However, 
the energies obtained with the Coulomb-Sturmian basis 
are also falling significantly more rapidly with increas- 
ing Amax- (In general, a higher starting energy for the 
convergence, at low N max , need not imply a lower rate 
of convergence.) Therefore, let us compare the exponen- 
tial fit parameters [see (62)] near the variational mini- 
mum. For the oscillator basis at Ml = 20 MeV, the con- 
vergence rate is c w 0.35, with an extrapolated ground 
state energy of —31.3 MeV. For the Coulomb-Sturmian 
basis at Ml — 30 MeV, the convergence rate is compa- 
rable, albeit marginally lower, at c w 0.29, with an ex- 
trapolated ground state energy which is also compara- 
ble, at —31.2 MeV. Interestingly, the extrapolations for 
the Coulomb-Sturmian basis have a qualitatively differ- 
ent dependence on Ml than those for the oscillator ba- 
sis. Rather than varying monotonically (increasing with 
increasing Ml), they have a minimum, at an Ml approx- 
imately equal to that of the variational minimum. 

The one significant numerical approximation which is 
entailed in setting up the Coulomb-Sturmian calcula- 
tions, as discussed in Sec. Ill C, is in the transformation 
of the two-body matrix elements of the nucleon-nucleon 
interaction from the oscillator basis to the Coulomb- 
Sturmian basis. The transformation is necessarily carried 
out in a truncated oscillator basis. It is therefore imper- 
ative to establish the numerical stability of the results 
with respect to the shell truncation N cu t in the sum over 
oscillator states. Calculations based on two-body matrix 
elements obtained with N cut — 9 (10 shells), N cut = 11 
(12 shells), and N cut — 13 (14 shells) are overlaid in 
Fig. 4(b), as well as in all subsequent plots of Coulomb- 
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FIG. 4: The 6 Li 1 + ground state energy (top) and 3 + excited state energy (bottom), calculated using the conventional harmonic 
oscillator basis (left) and the Coulomb-Sturmian basis (right). Calculated energies are plotted as a function of the basis Ml 
parameter, for iV max = 4, 6, 8, and 10 (successive curves, as labeled). For the Coulomb-Sturmian basis, calculations are shown 
variously for truncations iV cu t = 9 (dotted curves), iV cu t = 11 (dashed curves), and iV cu t = 13 (solid curves) in the change- 
of-basis transformation of two-body matrix elements. Exponentially extrapolated values (based on the iV C ut = 13 calculations 
in the case of the Coulomb-Sturmian basis) are indicated by crosses ( x ) . The best extrapolated values from the large-basis 
calculations of Ref. [35] are shown as horizontal dashed lines. 
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FIG. 5: The Li 1 + ground state energy, calculated using the 
Coulomb-Sturmian basis, but without making use of the sep- 
arable method of Sec. Ill D for the two-body matrix elements 
of the T rc i operator, for comparison with Fig. 4(b). That is, 
the entire Hamiltonian, including T re i, is transformed from 
the oscillator basis following the approach of Sec. IIIC. See 
the caption to Fig. 4 for further explanation of curves and 
symbols. 



Sturmian calculations. Calculations of the ground state 
energy for Ml > 20 MeV are highly stable with respect 
to this cutoff, in the present calculations. This range 
safely covers the variational minimum. However, the cal- 
culations are not stable with respect to this cutoff for 
Ml < 20 MeV, and higher cutoffs would therefore be re- 
quired for accurate results at these Ml values. The in- 
stability with respect to oscillator basis cutoff appears 
to increase with increasing iV max . Such a dependence is 
reasonable, since higher-A^ max calculations increasingly 
probe higher-n Coulomb-Sturmian single-particle basis 
functions, which in turn require a higher N cut for ac- 
curate expansion in an oscillator basis, as illustrated in 
Fig. 3. 

For the calculations shown in Fig. 4, the kinetic energy 
matrix elements have been calculated by the separable 
method of Sec. Ill D. It is interesting at this point to in- 
vestigate how essential it is to use the separable approach, 
rather than simply transforming the kinetic energy ma- 
trix elements from the oscillator basis. For comparison, 
we therefore repeat the calculations for the ground state 
energy in the Coulomb-Sturmian basis, but transforming 
the two-body matrix elements of the entire Hamiltonian 



from the oscillator basis, yielding the results shown in 
Fig. 5. It is seen that, without the separable calculation, 
the results are unstable with respect to N cut throughout 
the entire range of Ml values, including the vicinity of the 
variational minimum. Thus, the separable method plays 
a major role in obtaining numerically accurate calcula- 
tions. It would otherwise be necessary to start from oscil- 
lator two-body matrix elements in a significantly larger 
number of oscillator shells, possibly prohibitively so. 

Calculations for the energy of the 3 + first excited state 
for the oscillator basis are shown in Fig. 4(c) and for the 
Coulomb-Sturmian basis in Fig. 4(d). The results are 
very similar in nature to those for the ground state, so 
little additional discussion is required. The best extrap- 
olation from Ref. [35] places this state at 2.56(2) MeV 
excitation energy, corresponding to E w —28.93 MeV. 
For the oscillator basis at Ml — 20 MeV, the conver- 
gence rate is c w 0.34, with an extrapolated ground 
state energy of —28.8 MeV. For the Coulomb-Sturmian 
basis at Ml = 30 MeV, the convergence rate is again 
marginally lower, at c w 0.30, with an extrapolated en- 
ergy of —28.6 MeV, apparently erring on the high side 
relative to Ref. [35]. 

From these exploratory calculations for 6 Li, it would 
appear that convergence properties for energies with 
the Coulomb-Sturmian basis are comparable, i.e., not 
markedly inferior, to those of the oscillator basis, with 
some qualitative differences in the Ml dependence. We 
note that these exploratory results have not yet probed 
the variational freedoms available with the Coulomb- 
Sturmian basis, both in the choice of length parame- 
ters (Sec. IIIB) and in truncation schemes, as described 
above. The convergence rate alone does not provide con- 
clusive information on the robustness which can be ex- 
pected from large- _/V max extrapolation or on the best ex- 
trapolation procedure. Some questions regarding extrap- 
olation may be elucidated by extending the calculations 
to higher iV max . Furthermore, the rates of convergence 
of calculations with the oscillator and Coulomb-Sturmain 
bases will depend on the physical properties of the nu- 
cleus (and particular state) under consideration. For in- 
stance, the asymptotic properties of the single-particle 
basis may well play a larger role for halo nuclei or for 
states involving clusters with significant spatial separa- 
tion. 



C. Root-mean-square radius 

The root-mean-square (RMS) radius presents chal- 
lenges for convergence in NCCI calculations with the con- 
ventional oscillator basis [36]. Here we consider the in- 
trinsic, point- nucleo n RMS radius for the ground state, 
defined by \J (r^j) (see Appendix A), from which the 
center-of-mass contribution has been removed by con- 
struction. Evaluation of the expectation value (rf el ) in 
a many-body state requires that one first calculate the 
two-body matrix elements of the r^ el operator. These are 
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FIG. 6: The 6 Li 1 + ground state RMS radius, calculated using the conventional harmonic oscillator basis (left) and the 
Coulomb-Sturmian basis (right). Calculated energies are plotted as a function of the basis Ml parameter, for iV max = 4, 6, 
8, and 10 (successive curves, as labeled). For the Coulomb-Sturmian basis, calculations are shown variously for truncations 
N cu t — 9 (dotted curves), iV cu t = 11 (dashed curves), and iV cu t = 13 (solid curves) in the change-of-basis transformation of 
two-body matrix elements. Exponentially extrapolated values (based on the AT C ut = 13 calculations in the case of the Coulomb- 
Sturmian basis) are indicated by crosses ( x ) . The best estimated value from the large-basis calculations of Ref. [35] is shown 
as a horizontal dashed line. 



obtained for the Coulomb-Sturmian basis by the separa- 
ble method of Sec. Ill D. 

The oscillator basis results for the RMS radius in 
Fig. 6(a) are shown for the same range of calculations 
(JVmax = 4, 6, 8, and 10, with Ml values from 10 MeV to 
40 MeV) as for the energies in Sec. IV C. (The curves pro- 
ceed from greatest to least slope with increasing V max in 
the figure.) Exponential extrapolations to infinite N max 
are shown as well. The extrapolated values vary strongly 
with Ml and converge very slowly with V max . For in- 
stance, taking Ml at the variational minimum for the en- 
ergy, i.e., Ml w 20, the exponential convergence rate for 
the RMS radius with respect to iV max is only c w 0.024, 
and the extrapolated radius lies <~ 1 fm above the cal- 
culated values. Alternatively, the value at the crossover 
point of the curves obtained for different V max has also 
been proposed as an estimate of the full-space value [36] . 
This crossover occurs at Ml rts 12 MeV in the present cal- 
culations and lies in the vicinity of 2.2 fm. The best esti- 
mate from Ref. [35] , similarly obtained from the crossover 
point, for calculations with V max < 16, is ~ 2.3 fm, indi- 
cated by the dashed horizontal line in Fig. 6. 

Examining the calculations for the RMS radius using 
the Coulomb-Sturmian basis, as shown in Fig. 6(b), the 



gross features are similar. The crossover point for the 
curves obtained with iV max lies at Ml rts 20. The value 
of <~ 2.3 fm is consistent with the estimate of Ref. [35] 
and ~ 0.1 fm higher than the crossover for the curves 
obtained with the oscillator basis, for the same N max , in 
Fig. 6(a). Moreover, it is seen that exponential extrapola- 
tion may be a viable approach to estimating the full-space 
value for the radius. The extrapolated values obtained 
for Ml > 20 MeV, i.e., above the crossover point, are rea- 
sonably insensitive to Ml and are consistent with the best 
estimate from Ref. [35]. For instance, taking Ml w 30, 
i.e., at the variational minimum, the exponential con- 
vergence rate for the RMS radius is c ss 0.19, and the 
extrapolated radius is <~ 2.28 fm. Results are stable with 
respect to the shell cutoff in the transformation of matrix 
elements from the oscillator basis, for Ml > 20 MeV, as 
observed above for the energies. 

It would thus appear that the rate of convergence of the 
RMS radius obtained with the Coulomb-Sturmian basis 
is superior to that obtained with the conventional oscil- 
lator basis. However, further systematic investigation is 
required, especially into the stability of extrapolations 
with increasing V max , before general conclusions may be 
drawn. 
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FIG. 7: Expectation value of the number operator N c ul for 
center-of-mass oscillator quanta, as a function of oscillator en- 
ergy Ml. These calculations are for the 6 Li 1 + ground state, 
using the Coulomb- Sturmian basis, with Ml = 30 MeV. Cal- 
culations are shown for -/V max = 4, 6, 8, and 10 (successive 
curves, top to bottom), and for truncations iV C ut = 9 (dot- 
ted curves), iV cu t = 11 (dashed curves), and -/V cut = 13 (solid 
curves) in the change-of-basis transformation of two-body ma- 
trix elements. The analogous curve expected for a pure har- 
monic oscillator 0s function, with Ml — 17.5 MeV, is also 
shown for comparison (dotted curve, labeled). 



D. Center-of-mass dynamics 

We now focus on the dominant concern in using any 
basis other than the harmonic oscillator basis with A max 
truncation for nuclear many-body calculations, namely, 
incomplete separation of center-of-mass and intrinsic dy- 
namics. There are several aspects to consider: the natu- 
ral degree of separation arising in calculations using the 
Coulomb-Sturmian basis, the spurious state spectrum 
obtained in such calculations, and the extent to which 
a Lawson term can be used to influence spurious excita- 
tions. 

The problem of correcting for, or eliminating, spurious 
contributions for calculations with a general truncated 
basis is unresolved [37]. Nonetheless, it is still possi- 
ble that factorized wavefunctions might approximately 
be obtained in a truncated space. In the full space, fac- 
torization is obtained due to the separable Hamiltonian, 
albeit with degeneracies in the center-of-mass wavefunc- 
tions multiplying each intrinsic state (Sec. II B). There- 
fore, as larger truncated spaces are taken, approaching 
this full space, the structure of the eigenstates may be 
expected to converge towards such factorized structure. 
For instance, a high degree of factorization has been re- 
ported for coupled-cluster calculations in light nuclei [38] . 
Furthermore, introducing a Lawson term to the Hamilto- 
nian, as in (18), may serve to "purify" the eigenstates so 
that the motion more closely approximates 0s center-of- 



mass motion, as proposed by Gloeckner and Lawson [20] . 
This Lawson term also pushes eigenstates dominated by 
other center-of-mass excitations higher in the spectrum. 
However, caution must be exercised in such use of the 
Lawson term, since any improved (or, at least, more 
oscillator-like) description of center-of-mass motion may 
be obtained at the expense of the quality with which the 
intrinsic wavefunction is approximated [39] . 

A first indication of the degree of separation in the 
many-body eigenstate is provided by the expectation 
value of the A c m operator. This operator is defined, 
for an arbitrary center-of-mass harmonic oscillator en- 
ergy Ml, by 



N u = — ^ 



1 / P 2 



hfl V 2Am N 



+ 



Am N n 2 R 2 



—)■ 



(63) 



where the tilde serves to distinguish Ml from the basis 
Ml parameter. As noted by Hagen et al. [38], if sepa- 
ration occurs, as ^£>(rj;<Tj) = V'c.m.(R)V ; in(rij; <Tj), and if 
V'c.m.(R-) happens to be an oscillator 0s function, corre- 
sponding to some oscillator energy Ml, then the many- 
body wavefunction will have (N^ m ) = 0. Evaluation 
of the expectation value (N^ m ) requires that one first 
calculate the two-body matrix elements of P 2 and R 2 , 
and thence of Afp m . These are readily obtained for 
the Coulomb-Sturmian basis by the separable method 
of Sec. HID, so evaluation is straightforward. 

The expectation value (N^ m ) is shown as a function of 
Ml for the 6 Li 1 + ground state in Fig. 7, for the Coulomb- 
Sturmian basis calculation with basis Ml = 30 MeV and 
no Lawson term. The minimum value of (N^ m ) is ob- 
tained at Ml « 17.5 MeV, shifting gradually towards 
lower Ml, which corresponds to larger center-of-mass os- 
cillator length 6 c m . = [h/ {AniN^l)] 1 ! 2 , with increasing 
A max . (The location of the minimum also depends mod- 
estly upon the choice of basis Ml for the calculation, in- 
creasing with Ml.) The minimum value of (N^ m ) de- 
creases with increasing A max , but it appears to be con- 
verging towards a nonzero value of <~ 0.2. The fact that 
(N^ m ) values significantly less than unity are obtained 
in the calculations indicates that a 0s oscillator func- 
tion dominates the center-of-mass motion, and that an 
approximate separation of center-of-mass and intrinsic 
functions is spontaneously arising. However, the nonzero 
limit indicates that, as the full space is approached, the 
separated center-of-mass function is not strictly taking 
the form of a 0s oscillator function. For comparison, the 
dependence of (N^ m ) on Ml which would be obtained for 
a pure oscillator 0s function with Ml = 17.5 MeV, given 
by (A^m.) = + H/O - 2), is also shown in Fig. 7. 

In interpreting these results, it must be stressed that 
calculating (N^ m ) for an eigenstate provides only a lower 
limit on the degree of factorization. That is, a nonzero 
(A^ m ) does not preclude factorization but can simply 
indicate that the factorized center-of-mass wavefunction 
is not of 0s oscillator type. Extracting the true degree of 
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FIG. 8: Level spectrum for 6 Li, including spurious states, calculated using the conventional harmonic oscillator basis. Calcula- 
tions (left to right) are for iVmax = 4, 6, and 8, and then shown again for iVmax = 8 with addition of a Lawson term, of sufficient 
strength to shift the spurious states above the energy range displayed in this plot. For each level, the angular momentum J 
is indicated at left, and (iV c . m .) is indicated at right. For degenerate multiplets of spurious states, the thickness of the line is 
proportional to the number of states. These calculation are for Ml — 20 MeV. 



factorization is more challenging. To do so likely requires 
some form of explicit transformation to center-of-mass 
and relative coordinates. For instance, an expansion ip = 

J2i Siipch.ipin mav then be obtained through a singular 
value decomposition as proposed in Ref. [38]. 

Since factorization arises in the full space, the effects 
of convergence of the ccntcr-of-mass dynamics were al- 
ready implicitly included in the extrapolations to the full- 
space values of the observables of interest, as explored in 
Sees. IV B and IV C. However, for this extrapolation to 
be possible, it is necessary that states involving spurious 
excitations of the center-of-mass function can be disen- 
tangled and removed from the low-lying spectrum. This 
becomes an increasing concern with increasing /V max , as 
we shall now see from examining the spurious state spec- 
trum. 

It is helpful to first consider the eigenvalue spectrum, 
including spurious states, obtained in calculations with 
an iV max -truncated oscillator basis. This is illustrated for 
6 Li in Fig. 8, for N max = 4, 6, and 8, with Ml = 20 MeV. 
The eigenvalue of N c m is indicated to the right of each 
level in this figure. For instance, at /V max = 4, where the 
two lowest states in the spectrum are shown, these states 
have 7V C m = and are thus nonspurious, corresponding 
to the intrinsic 1 + ground state and 3 + first excited state, 
with 0s ccntcr-of-mass motion. 



Let us examine the evolution of the spectrum in Fig. 8 
with increasing /V max , bearing in mind the direct sum 
structure (21) of the 7V max -truncated space. Moving 
from iV max = 4 to N max = 6, the two additional oscil- 
lator quanta introduced to the system may go towards 
converging the intrinsic states. This yields the new 1+ 
ground state and 3 + excited state at lower energies in 
the -/V max = 6 calculation (in the W° m . ® Hf n subspace). 
Alternatively, the two additional quanta may go into 
center-of-mass excitation, yielding spurious states (in the 
%c.m. ® Win subspace). Since the center-of-mass excita- 
tion gives no contribution to the energy, under the in- 
trinsic Hamiltonian we are using, the resulting spurious 
states are degenerate with the nonspurious states ob- 
tained at jV max = 4 (in the "H^.m. ® subspace) . Then, 
moving to iV max = 8, new N cm = 2 spurious states ap- 
pear degenerate with the /V max — 6 nonspurious states, 
new 7V C m =4 spurious states appear degenerate with 
the iV c m = 2 spurious states from iV max = 6, etc. 

To ascertain the angular momenta expected for the 
spurious states, we note that angular momentum cigen- 
states of the full eigcnproblem are obtained from those of 
the intrinsic eigcnproblem via angular momentum cou- 
pling as = [i/>c'm m ) x Vi ( n m) ] (J) - Thus, the angu- 
lar momenta expected for the spurious states follow by 
the triangle inequality for addition of the center-of-mass 
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FIG. 9: Level spectrum for 6 Li, calculated using the Coulomb-Sturmian basis. Calculations (left to right) are for iV max = 4, 
6, and 8, and then again for -/V ma x = 8 with addition of a Lawson term aN^m. of strength a = 2 MeV. Energies corrected by 
— a(iV^m.) are shown at far right. For each level, the angular momentum J is indicated at left, and (-/V c n m .) is indicated at 
right, as a measure of the number of center-of-mass oscillator quanta. Approximately degenerate multiplets of spurious states 
are marked by brackets and connected to the state at lower JV max to which they are approximately related by coupling to two 
center-of-mass quanta. The dashed lines trace the change in level energy induced by the Lawson term. For JV max = 8, an arrow 
connects the two J — 2 levels which may be described (see text) as admixtures of a nonspurious and spurious level. These 
calculation are for HQ = 20 MeV, with 7V cut = 13. The quantity (A r ,P m ) is evaluated for M7 = 20 MeV, and the Lawson term is 
defined for W. L = 20 MeV as well. 



angular momentum l c , m , and intrinsic angular momen- 
tum Ji n . Recall that the three-dimensional oscillator 
spectrum contains angular momenta I = for N = 0, 
I = 1 for N = 1, I = (0,2) for N = 2, I = (1,3) for 
N = 3, I = (0,2,4) for N = 4, etc. Spurious states with 
N c m = 1 lie in the opposite-parity space and therefore 
do not appear in Fig. 8. 11 However, for N c m — 2, cou- 



pling Z c . m . = and 2 to the J; n = 1 intrinsic ground 
state yields a spurious-state multiplet with angular mo- 
menta (3,2,1,1), as seen in Fig. 8. Similarly, coupling 
these values of Z c . m . to the J; n = 3 intrinsic excited state 
yields a spurious-state multiplet with angular momenta 
(5,4,3,3,2,1). The N c m = 4 spurious multiplets seen 



the 'Kc.m. ® %f n su b s P ace ! do indeed appear in the even-parity 
Odd spurious exitations of odd-parity intrinsic states, e.g., in spectrum, but in 6 Li these are at higher energy. 
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FIG. 10: The 6 Li 3+ excited state energy, calculated using the 
Coulomb-Sturmian basis, as in Fig. 4(d), but now including 
a Lawson term, with strength a — 2 MeV and Lawson term 
oscillator energy KIl chosen equal to the basis Ml. Energies 
are corrected by subtracting a(JV c Q m.). See the caption to 
Fig. 4 for further explanation of curves and symbols. 



at A max = 8 in Fig. 8 are obtained similarly by coupling 
^c.m. — 0j 2, and 4 to the intrinsic state. 

It is apparent from Fig. 8 that a high level of contami- 
nation of the low- lying spectrum with spurious states will 
arise with increasing A max , as the difference in energy be- 
tween intrinsic ground states in successive -/V max spaces 
decreases. Already several spurious states arise below 
the first excited state at A max = 8 in this example. This 
would become a prohibitive problem for large A max , as 
suggested in Sec. II C, since the Lanczos diagonalization 
in the many-body problem must converge the spurious 
states along with the nonspurious states. However, for 
the A max -truncated oscillator b noted in Sec. II C, 

inclusion of the Lawson term in the Hamiltonian pushes 
the spurious solutions to higher energy, without affecting 
the nonspurious states, obviating the problem [Fig. 8 (far 
right)]. 

With this understanding of the spurious state spec- 
trum for the oscillator basis, we now have a baseline 
for interpreting the eigenvalue spectrum obtained with 
a Coulomb-Sturmian basis, shown in Fig. 9 for basis 
Ml = 20 MeV. It is seen that the same multiplets of spu- 
rious states (marked with brackets in the figure) arise as 
in the calculation based on the oscillator basis, but now 
the degeneracies — with the nonspurious state at lower 



AT max and between the members of the multiplct itself — 
are only approximate. Since we are using the intrinsic 
Hamiltonian, these energy differences do not arise from 
any direct contribution of the center-of-mass dynamics 
to the energy. Rather, to the extent that factorization 
occurs, these differences arise from variation in the level 
of convergence of the intrinsic wavefunction associated 
with the center-of-mass wavefunction. Alternatively, to 
the extent that factorization is imperfect, these differ- 
ences can arise from admixtures of contributions involv- 
ing different ccntcr-of-mass and intrinsic excitations. 

Although these states in Fig. 9 are not eigenstates of 
N^ m , we can still calculate an average number of center- 
of-mass oscillator quanta as (N^ m ) . This expectation 
is indicated to the right of each level in Fig. 9, where 
Ml has been chosen simply equal to the basis Ml 1 i.e., 
Ml = 20 MeV. It is seen that the (N^ m ) values clearly 
reflect the identification of the states as spurious or non- 
spurious according to the energy spectrum noted above. 
The nonspurious states share a similar range of values 
for (N^ m ) — at A max = 8, nearly identical for the 
ground state and first excited state (~ 0.24) and some- 
what higher (~ 0.3-0.4) for some of the higher states. 
The states analogous to the iV c m = 2 spurious states of 
the oscillator-basis calculation, in contrast, have (N^ m ) 
values which cluster closely around 2.4. 

Note the two 2 + states at about —21 MeV in the 
A max = 8 calculation of Fig. 9. With exact factorization, 
one of these would be nonspurious and the other spu- 
rious. However, the {N^ m ) values for these two states 
(~ 0.96 and ~ 1.74) indicate that the spurious and non- 
spurious states are strongly mixed. This mixing provide 
an illustration of the challenge associated with contam- 
ination of the low-lying spectrum with spurious states. 
As the density of spurious states increases with N max , 
the close proximity of spurious and nonspurious states 
may be expected to lead to extensive mixing and conse- 
quently a breakdown of center-of-mass factorization for 
even the lowest-lying states. Therefore, it is even more 
important that the spurious states be eliminated from the 
low-lying spectrum than it is for conventional oscillator 
basis calculations. 

With this in mind, we explore the efficacy of the Law- 
son term when used with the Coulomb-Sturmian basis. 
At right in Fig. 9, the effect of introducing a Lawson term 
to the Hamiltonian is shown. For simplicity in 
this illustration, we choose MIl — 20 MeV, correspond- 
ing to the basis Ml and the Ml for the (N^ m ) values 
indicated. This choice is arbitrary, 12 and another value, 



When used with the oscillator basis, the center-of-mass oscillator 
energy HQ^ in the Lawson operator is generally chosen equal to 
the basis HQ, to preserve factorization. However, when used with 
a general, non-oscillator basis, there is no such requisite pairing, 
and HQ l may be chosen freely, so as to obtain the most effective 
removal of spurious dynamics. 
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such as that at which ) is minimized for the ground 

state, might well be profitably used. A Lawson strength 
a = 2 MeV has been adopted, as sufficiently large to ex- 
punge spurious states from the lowest few MeV of the 
spectrum, but not so large as to place undue weight on 
coercing the center-of-mass wavefunction into a pure Os 
oscillator state, at the possible expense of compromising 
convergence of the intrinsic state. The change of energies 
with introduction of the Lawson term is traced by dashed 
lines in Fig. 9. Notice that the mixing of the nonspurious 
and spurious 2+ states discussed above (now at energies 
of about —20 MeV and —16 MeV, respectively) has been 
eliminated. 

The Lawson term is also seen, from the expectation 
values indicated in Fig. 9, to reduce (N^ ) for these 
states. It is not yet clear how much of this change reflects 
improvement of the center-of-mass factorization and how 
much simply relects modification of an already-factorized 
center-of-mass function towards oscillator 0s form. 

Since even the nonspurious states have nonzero val- 
ues for (N^ ), their energies are raised by introduction 
of the Lawson term, by ~ a(N^ ). This contribution 
is not expected to vanish in the large -/V max limit, since 
(A^iPm ) has already been seen not to approach zero. To 
recover the eigenvalue of the intrinsic Hamiltonian on 
the intrinsic wavefunction, to the extent that good fac- 
torization is obtained, we must correct the calculated 
energy for the contribution of the Lawson term acting 
on the center-of-mass function, by subtracting a(N^ ) 
back off. The energies obtained after this correction, for 
the nonspurious states, are shown at far right in Fig. 9. 
After correction, the original values for the energies, as 
obtained before introduction of the Lawson term, are al- 
most (but not quite) recovered. The corrected energies 
are still marginally higher, likely reflecting the compro- 
mise in convergence of the intrinsic state incurred by the 
Lawson term, and this discrepancy increases with the 
Lawson strength a used in the calculation. 

The Lawson term thus appears to be a credible means 
of eliminating spurious states from the low-lying spec- 
trum, in calculations with the Coulomb-Sturmian basis. 
The essential question, if the Lawson term is to be used in 
practice, is whether or not the Lawson term has any sig- 
nificant adverse impact on convergence properties. Tak- 
ing the energy of the first excited state as an example, 
we repeat the calculations of Fig. 4(d), but now including 
a Lawson term of strength a = 2 MeV, resulting in the 
energies in Fig. 10. The a(N^ ) correction to the en- 
ergies, described above, has been included. The results 
are virtually indistinguishable from those of Fig. 4(d). 
For comparison with the discussion in Sec. IV B, we note 
that the convergence rate at the variational minimum 
{Ml = 30 MeV) is still c w 0.30, and the extrapolated 
energy is still approximately —28.6 MeV. 



V. CONCLUSION 

Although the conventional oscillator basis has definite 
advantages for ab initio nuclear many-body calculations 
with the NCCI approach, namely, the potential for ex- 
act center-of-mass factorization of eigenstates and the 
simplicity of the Moshinsky transformation for Hamil- 
tonian matrix elements, it also presents the disadvantage 
of nonphysical Gaussian asymptotics at large distances, 
i.e., the oscillator wavefunctions satisfy the wrong bound- 
ary conditions at infinity for use with bound states of 
nuclei. The Coulomb-Sturmian functions retain the ad- 
vantages of forming a complete, discrete set of square- 
intcgrablc functions while also exhibiting realistic ex- 
ponential asymptotics. We have seen that the techni- 
cal and physical challenges of carrying out NCCI cal- 
culations with a Coulomb-Sturmian basis are tractable. 
To briefly summarize the computational framework, the 
many-body calculation has the standard structure for 
an nlj single-particle basis, the interaction matrix ele- 
ments are transformed from the harmonic-oscillator ba- 
sis, and relative kinetic energy matrix elements are cal- 
culated separably. In the initial exploratory calculations 
considered here, it is found that the convergence rates 
for energies are competitive with those obtained with an 
oscillator basis, the convergence rate for the RMS ra- 
dius is superior, and spurious center-of-mass excitations 
can be successfully managed. Many of the considera- 
tions addressed in this work could be relevant to NCCI 
calculations with other possible radial bases as well, e.g., 
transformed harmonic oscillator bases [40] . 

The importance of the asymptotic properties of the 
basis functions may be expected to vary depending upon 
the physical properties of the nucleus and state under 
consideration. A basis such as the Coulomb-Sturmian 
basis might well be particularly appropriate to halo nu- 
clei, where the mismatch with the oscillator functions at 
large distances is particularly severe. Another case of in- 
terest would be states involving clusters with significant 
spatial separation. The importance of reproducing the 
large-r properties of the nuclear eigenstates may also be 
expected to depend upon the observable under consider- 
ation, depending upon how heavily large-r contributions 
are weighted by that observable. Thus, e.g., the differ- 
ence between Gaussian and exponential asymptotics may 
be expected to be more important for the RMS radius or 
E2 observables than for Ml observables. Asymptotic 
properties also play a significant role in scattering prob- 
lems. The extent to which a Coulomb-Sturmian basis 
may be successfully used in ab initio scattering calcula- 
tions, e.g., through a generalization of the no-core shell 
model resonating group method (NCSM/RGM) [41], will 
depend critically upon the details of the center-of-mass 
factorization properties. 

To more fully ascertain the relative advantages or dis- 
advantages of the Coulomb-Sturmian basis for NCCI cal- 
culations, extensive and systematic calculations are re- 
quired, both into the convergence properties of the basis 
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and the robustness of extrapolations. Most obviously, 
these need to be carried to high -/V max , for a variety 
of nuclei and interactions. However, there is also con- 
siderable room for optimization within the method it- 
self, which must be explored. The prescription for the 
l-dependence of the length parameter within the single- 
particle basis (Sec. Ill B) and the many-body truncation 
scheme (Sec. IV A), in which the present oscillator-like 
N = 2n + I "energy" weighting is dictated purely by 
convenience, are notable areas of possible improvement. 
Although the two-body JISP16 interaction was used in 
the illustrative calculations, the transformation proce- 
dure (Sec. IIIC) carries over readily to three-body in- 
teractions, so convergence properties with, e.g., chiral 
effective field theory interactions with similarity renor- 
malization group evolution [42], can be investigated. 
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Appendix A: Center-of-mass decomposition of r 2 
and k 2 



The one-body operators r 2 = J2i r i an d k 2 = J^fe^, 
for the A-body system, may be decomposed into separate 
parts depending only upon the center-of-mass coordinate 
(or momentum) and on the relative coordinates (or mo- 
menta), respectively. The decompositions of the kinetic 
energy operator T and noninteracting harmonic oscillator 
potential U into center-of-mass and relative parts follow 
immediately. In this appendix, we summarize the rela- 
tions among relative and center-of-mass operators, both 
for reference in the present discussion and to establish a 
uniform notation for the description of coordinate-space 
and momentum-space matrix elements in Sec. Ill D. 

Recall that the center-of-mass coordinate and momen- 
tum vectors are 



R 



(Al) 



In the following, we let p 4 = hki, P = KK, etc. 

First, consider the one-body r 2 operator, defined by 



= £ 



r 2 



(A2) 



Comparing the sum on the right hand side of (A2) with 
those in the operators 13 

i i ij 

A 2 r 2 el = I £'(r* - r,) 2 = (A - 1) £ r 2 ^'r, • Tj . 



demonstrates that 



2 _ 4 2 d 2 i \1 2 

r rcl- 



Ar z = A Z R Z + A 



(A3) 



(A4) 



Multiplying by (mf2 2 )/(2A) gives the decomposition of 
the harmonic oscillator potential energy operator U n into 
center-of-mass and relative contributions, U n = U^ m + 
U^, where 



U 



2 D 2\ pfi 

mQ 2 



(A5) 



2A 



{At 2 ). 



The quantity r 2 el has the geometric significance that it 
is the mean square radius relative to the center of mass, 
i.e., 



(A6) 



The square root of the expectation value of this operator, 
(r 2 ^) 1 / 2 , is the point-nucleon RMS radius. 

Similarly, consider the one-body k 2 operator, defined 

by 



k 2 = J2kl 



(A7) 



Comparison with the sums in 



K 2 = £ k 



— k 2 + ^ ki ■ kj , 



k 2 cl = \ ]T'(k,> - k,) 2 =(A-l)Y^ k 2 - • k 3 , 

ij i ij 

(A8) 



13 We include the factors of A 2 on the left hand side of (A3) as com- 
pensation for the factor of l/A appearing in the definition (Al) of 
R., so as to simplify the right hand side. In particular, this main- 
tains the parallel with the decomposition of momentum space 
operators in (A8) 
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demonstrates that 



Ak 2 = K 2 + k 2 cl . 



(A9) 



Multiplying by ft 2 /(2Aitin) gives us the decomposition 
of the kinetic energy operator T into center-of-mass and 
relative contributions, T — T c . m . + T rc i, where 



(ft 2 * 2 ) Trel= (^4i) ! 
2 Am at ' rc 2AniN ' 
_ (Afi 2 fc 2 ) 



(A10) 
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Appendix B: Zeros of generalized Laguerre and 
Jacobi polynomials 

Numerically robust evaluation of the radial inte- 
grals which arise in evaluation of the overlaps be- 
tween harmonic oscillator and Coulomb-Sturmian bases 
(Sec. Ill C) and the radial matrix elements for the 
Coulomb-Sturmian basis (Sec. HID) requires accurate 
knowledge of the zeros of the integrands in (40)-(41), 
(49)-(50), and (58)-(59), thus of generalized Laguerre 

polynomials L"(x) and Jacobi polynomials J n a ^\x). Al- 
though, in principle, generic numerical rootfinding algo- 
rithms may be used, it is preferable to determine the 
zeros according to a more reliable and efficient approach 
specific to orthogonal polynomials, such as the Golub- 
Welch algorithm [43]. This method requires recurrence 
coefficients for the relevant monic polynomials, i.e., such 
that the highest-order coefficient is unity, as summarized 
in this appendix. 

The Golub- Welch algorithm is specifically formulated 
for monic polynomials. Consider a family of polynomials 
Pn(x) = J2m=a c m.x m (n = 0, 1, ...), orthogonal un- 
der weight function w(x) on the interval [a,b], and with 
c„ = 1. Suppose these polynomials satisfy the recurrence 
relation 

Pn+l{x) + (B n - X)p n (x) + A n p n _ 1 (x) = 0, (Bl) 

characterized by recurrence coefficients A n and B n . 
Then, to find the zeros p n via the Golub- Welch algo- 
rithm [43], one must construct the corresponding Jacobi 
matrix J. This is the n x n tridiagonal matrix consist- 
ing of entries Jis — on the main diagonal and 
Ji-i t i — Ji,j-i = (Ai-i) 1 / 2 on the adjacent diagonals. 
As a tridiagonal matrix, J is easily diagonalized. The 
eigenvalues Xi, for i = 1, 2, . . ., n, are then the zeros of 

Pn- 

The generalized Laguerre polynomials L™ are not 
monic, having c„ = (— )"n! [28]. We must therefore in- 
stead consider the monic generalized Laguerre polyno- 
mials L°, defined by L°(x) = (-) n n\L%(x) [44]. These 
satisfy a recurrence relation of the form (Bl), with recur- 
rence coefficients 



A n = n(n + a) 
B n = 2n + a + 1. 



The Jacobi polynomials P n a '^ are likewise not monic, 
having c„ = 2-"( 2n+ n Q+/3 ) [28]. We must therefore in- 
stead consider the monic Jacobi polynomials P n a ' l3 \ de- 
fined by Pi a ' 0) (x) = 2 n ( 2n+ ^ +!i y 1 P n a ' fi \x) [44]. These 
satisfy a recurrence relation of the form (Bl), now with 



A n 



B„ — 



4n(n + a)(n + f3)(n + a + f3) 



(2n + a + p) 2 {2n + a + /3 + l)(2n + a + /3 - 1) 

!3 2 -a 2 

(2n + a + /3)(2n + a + /3 + 2) ' 



(B2) 



(B3) 

The Golub- Welch algorithm also yields the weights 
Wi appearing in the n-point Gaussian integration 
formula associated with this family of polynomials, 

j^f{x)w{x)dx w Y^i=\ f( x i) w iy which are obtained 
from the eigenvectors of J, as detailed in Ref. [43]. The 
zeros and weights appearing in n-point Gauss-Legendre 
quadrature formulas used in evaluating the radial inte- 
grals of Sec. Ill are widely tabulated [28]. However, it is 
convenient to note that the Jacobi matrix required in ob- 
taining these may also be obtained using the recurrence 
coefficients of (B3). It is necessary to consider the monic 
Legendre polynomials P n (x) = [2™ (n!) 2 / {2n)\]P n (x) [44], 
which constitute a special case of the monic Jacobi poly- 
nomials, P n (x) = P n '°\x), described by (B3) with 
a = P = 0. 



Appendix C: Two-body states 

In this appendix, the notation is established for the 
antisymmetrized (AS) and normalized antisymmetrized 
(NAS) two-particle states, with angular momentum cou- 
pling. These definitions are required for discussion of 
like-particle two-body matrix elements in Sec. III. 

We first define angular momentum coupled states 

\ab;J) = ^ (j a m a j b m b \JM)\am a )\bm b ). (CI) 



for two distinguishable particles, that is, particle 1 is 
in orbital a, and particle 2 is in orbital b. We denote 
such distinguishable-particle states by using parentheses 
rather than angle brackets, following the conventions of 
Ref. [29]. Such states may be used directly for the case 
of one proton and one neutron, i.e., 

\ab;J) pn = ^ (j a m a j b m b \JM)\am a ) p \bm b ) n . (C2) 



However, for two like fermions, antisymmetrized states 
are then obtained as 

\ab; JM) AS = (4x C t)^|) 

= -L[\ab;JM) - (-) J -^\ba;JM)]. 

(C3) 
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These states have the basic symmetry property 

\ab; JM) As = -(-) J - j *- j *\ba; JM) as- (C4) 

Therefore, if the orbitals a and b are identical, only states 
with J even may be obtained. The states defined in (C3) 
are antisymmetrizcd but not strictly normalized, in that 
a further factor of l/y/2 is required for normalization in 
the special case in which both particles occupy the same 
orbital. Strict normalization, even in this special case, is 
obtained by taking normalized antisymmetrized (NAS) 
states 

\ab; JM) NAS = (1 + 6 ab )-^ 2 \ab; JM) AS - (C5) 

Two-body matrix elements may be represented in ei- 
ther the AS scheme or NAS scheme, with the relation 

(cd; J\V\ab; J) N as 
= (1 + M- 1/2 (l + 5 ab )- 1/2 (cd- J\V\ab; J) AS , (C6) 

shown here for matrix elements of a scalar operator V 
within a single J-space. Both schemes are in common 
use for representing interaction matrix elements. The 
AS scheme may yield simpler expressions than the NAS 
scheme, e.g., as seen comparing the change of basis rela- 
tion (45) with (46). 



and (cd; J\ Vki-k 2 \ a b; J)o, by which we denote matrix ele- 
ments evaluated for &ho = 1. In this appendix, we give 
explicit expressions for matrix elements of physically rel- 
evant operators, for a given basis Ml, in terms of these 
reference matrix elements and dimensional scale factors. 

Two-body matrix elements of the relative kinetic en- 
ergy are given by 



[Ml\ 

(cd;J\T Icl \ab;J) = ( — J (cd; J\V k 2 \ab; J) 

-2(^j(cd;J\V kl . k2 \ab;J) (Dl) 
and those of the r 2 el observable by 

(cd; J\r 2 cl \ab; J) = (^^j {cd; J\V r ,\ab; J) 
bio 



-2[^r)(cd;J\V ri . r2 \ab;J) Q . (D2) 

For present purposes, it is most convenient to reexpress 
&ho in terms of Ml using combinations of physical con- 
stants chosen so as to only involve energy and length 
units, as 



(he) 



[(m N c 2 )(hn)] 1 / 2 ' 



(D3) 



Appendix D: Rescaling of separable matrix elements 

For the separable calculation of matrix elements 
described in Sec. HID, the relations (61) pro- 
vide A-dependent expressions for the two-body ma- 
trix elements of R 2 , r 2 cl , K 2 , and k 2 cl in terms 
of the A-indepcndcnt two-body matrix elements 
(cd; J\V r 2 \ab; J), (cd; J\V Tl . T2 \ab; J), (cd; J\V k 2 \ab; J), 
and (cd; J\V kl .^ 2 \ab; J). These matrix elements still de- 
pend on the length parameter chosen for the basis. How- 
ever, the operators V r 2 and V ri . T2 are homogeneous of 
order 2 in the coordinates, i.e., their matrix elements 
scale with the length parameter as b 2 , and the opera- 
tors Vfe2 and Pk r k 2 are homogeneous of order —2, i.e., 
their matrix elements scale as b~ 2 . Recall that, under 
the prescription of Sec. IIIB, the length parameters bi 
appearing in all Coulomb-Sturmian functions are pro- 
portional to a common length parameter &ho (this com- 
mon proportionality is a general property to be expected 
of any prescription for the bi). Therefore, these matrix 
elements of V r i, V Tl . r2 , V k 2, and Vki-k 2 need only be cal- 
culated once, at some particular reference value for the 
length scale, and then may be transformed to the actual 
length scale, or Ml value, of the many-body calculation 
by simple multiplication. For evaluation of the radial in- 
tegrals appearing in (49), (50), (58), and (59), it is nat- 
ural to adopt a dimensionless reference scale &ho = 1- 
Thus it is only necessary to evaluate matrix elements 
(cd; J\V r 2\ab; J) , (cd; J\V ri . r2 \ab; J) (cd; J\V k 2\ab; J) , 



where m N c 2 w 938.92 MeV and he w 197.327MeVfm. 

In the investigation of center-of-mass separation and in 
the Lawson term as applied to NCCI calculations with 
the Coulomb-Sturmian basis (Sec. IV D), we consider the 
center-of-mass oscillator number operator ATp m of (63), 
involving an arbitrary oscillator energy Ml, in general 
different from the basis Ml. This operator has two-body 
matrix elements 

(cd;J\N ( 2 m \ab;J) = 
1 {Ml) 



2A (Ml) 

i (hCi) 



2A (Ml) 



A l (cd; J\V k 2\ab; J) +2(cd; J\V kl . k2 \ab; J) 

1 - (cd; J\V r 2\ab; J) +2(cd; J\V Tl . T2 \ab; J) 
--^—^(cd;J\t 2b \ab;J), (D4) 



where 1 2 & is the identity operator on the two-body space. 
If one is evaluating (N^ m ) for several values of Ml, as 
in Fig. 7, it suffices to calculate the expectation values 
of just the two operators P 2 and R 2 for the many-body 
state, since these two numerical values may then be com- 
bined arithmetically by (63) to deduce (N^ m ) for any 
value of Ml. More simply, in terms of the expectation 
values of the dimensionless operators Kq = P 2 /(m^Ml) 
and A 2 Rq = (m^VlA 2 /h)R 2 , corresponding to the two- 
body matrix elements appearing in brackets, respectively, 
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in (D4), we have 



i m 2 i («j) 



(A 2 i?g)--. (D5) 



The number operator N n for the one-body harmonic 
oscillator Hamiltonian H n , though not used in the 
present work, can also be of interest. For instance, if a 
many-body state has been obtained from an NCCI calcu- 
lation using the Coulomb-Sturmian basis, (N Q ) provides 
an estimate of the number of quanta which would be re- 
quired to represent this same state in a space spanned by 
a conventional harmonic-oscillator basis of oscillator en- 



ergy Ml, or a Hamiltonian term proportional to N may 
be used for calculations involving an external harmonic 
oscillator trapping field. The two-body matrix elements 
are given by 



1 (nn) 

2(A - 1) (ftft) 
1 (fth) 



(cd;J\N%b;J) = „,/ )^y (cd;J\V k 2\ab; J) 

J- J (nil) 



2(A-1) (Ml) 



(cd; J\V r 2 \ab; J)o 

(cd;J\t 2b \ab;J). (D6) 
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